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Under  AFOSR  Contract  F49620-87-K-0003  which  was  initially  funded  November  1, 


1986  and  completed  on  September  30,  1990,  the  problem  of  control  augmented  structural 
optimization  of  aeroelastically  tailored  fiber  composite  wings  was  addressed  in  a  series  of 
comprehensive  studies.  This  research  culminated  in  the  first  truly  integrated,  practical 
computer  program  capable  of  treating  this  multidisciplinary  synthesis  problem  by 
simultaneously  changing  structural,  aerodynamic  and  control  type  design  variables  for 
practical  aircraft  configurations  such  as  the  F-16  fighter. 


The  main  line  of  this  research  program  has  been  described  in  a  substantial  number 
of  publications  [1-9].  A  brief  description  and  summary  of  the  contributions  contained  in 
these  publications  is  provided  below.  A  more  detailed  description  of  this  research  effort  is 
summarized  in  Refs.  3,  4,  6  and  8  which  are  attached  as  Appendix  A  of  this  report. 


Reference  1  is  a  detailed  report  describing  the  structural  model  used  in  this  research 
activity  and  its  implementation  in  a  computer  code.  Using  this  structural  model  a  complete 
airplane  configuration  can  be  efficiently  modeled  as  an  assembly  of  flexible  lifting  surfaces. 
Each  lifting  surface  is  modeled  as  an  equivalent  plate  whose  stiffness  is  controlled  by 
contributions  from  thin  cover  skins  (fiber  composite  laminates)  and  the  internal  structure 
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(spar  and  rib  caps).  Wing  sections  are  connected  to  each  other  via  stiff  springs 
(representing  hinge  stiffnesses  at  attach  points)  and  flexible  springs  (representing  the 
stiffness  of  actuators  and  their  backup  structure).  Each  wing  section  can  include  several 
trapezoidal  parts.  Concentrated  masses  are  used  to  model  nonstructural  items  and  balance 
masses. 

Structural  topology,  shape  and  material  properties  are  preassigned;  however,  skin 
layer  fiber  orientations  are  treated  as  design  variables.  Skin  thicknesses  are  also  indirectly 
represented  as  design  variables.  Concentrated  masses  and  springs  at  preassigned  locations 
can  also  be  treated  as  design  variables.  It  was  demonstrated  that  this  structural  model  is 
capable  of  capturing  the  important  modes  of  an  F-16  fighter  wing,  with  accuracy  comparable 
to  a  detailed  finite  element  model,  and  at  the  fraction  of  the  computational  cost.  This 
computational  efficiency  is  a  key  ingredient  in  the  success  of  the  subsequent  optimization 
studies,  based  upon  this  model. 

Reference  3,  which  is  an  expanded  and  improved  version  of  Ref.  2,  presents  the 
theoretical  basis  for  the  synthesis  of  an  actively  controlled  composite  wing  as  a 
multidisciplinary  optimization  problem.  A  unique  integration  of  analysis  techniques 
spanning  the  disciplines  of  structures,  aerodynamics,  and  controls  is  described.  A  rich 
variety  of  behavior  constraints  can  be  treated  including  stress,  displacement,  control  surface 
travel  and  hinge  movement,  natural  frequency,  aeroservoelastic  stability,  gust  response,  and 
handling  quality  constraints,  as  well  as  performance  measures  in  terms  of  drag/lift 
coefficients,  drag  polar  shape,  required  load  factor  or  roll  rate,  and  wing  mass.  The  design 
space  includes  a  simultaneous  treatment  of  structural,  aerodynamic,  and  control  system 
design  variables.  The  basis  for  multidisciplinary  wing  optimization  is  prepared  by 
formulating  the  analysis  capability  together  with  the  related  behavior  sensitivity  analysis. 
Applicability  of  approximation  concepts  to  this  particular  multidisciplinary  optimization 
problem  is  examined  by  studying  typical  aeroservoelastic  stability,  gust  response  and 
performance-related  constraints.  The  computational  efficiency  of  the  combined  analysis  and 
sensitivity,  as  well  as  the  quality  of  key  behavior  constraint  approximations,  was  shown  to 
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be  excellent.  Thus  the  single-level  optimization  of  composite,  actively  controlled  practical 
wings, ^described  in  the  subsequent  references  could  be  carried  out. 

The  integrated  multidisciplinary  synthesis  capability  (described  in  Ref.  3)  was  used 
in  Ref.  4  in  a  series  of  innovative  exploratory  design  studies  in  which  constraints  from 
several  disciplines  are  taken  into  account  simultaneously  and  the  design  space  is  opened  up 
to  include  structural,  control  system,  and:  aerodynamic  design  variables.  The  effectiveness 
and  the  efficiency  of  the  new  capability  are  studied  using  a  mathematical  model  of  a 
remotely  piloted  vehicle  (RPV).  The  emphasis  in  these  studies  was  on  active- 
control /structure  interaction  in  the  RPV  wing  design  problem.  These  studies  have 
demonstrated  the  successful  adaptation  of  NLP/AC  (nonlinear  programming/approximation 
concepts)  techniques  to  problems  with  important  new  types  of  constraints  (aeroservoelastic 
stability,  gust  response)  exhibiting  greater  complexity  than  those  previously  treated  in  pure 
structural  synthesis  problems. 

In  Ref.  6  the  effectiveness  and  efficiency  this  integrated  aeroservoelastic  optimization 
capability  is  displayed  by  applying  it  to  an  RPV  as  well  as  more  complex  F-16  and  X-29  type 
airplane  models.  Simplified  handling  quality  constraints  are  added  to  the  set  of  design 
requirements.  The  performance  of  several  complex  eigenvalue  approximations  was  also 
examined.  Effects  of  control  law  structure  on  the  weight  and  robustness  of  the  resulting 
aeroservoelastic  design  provide  new  insights  into  the  complex  multidisciplinary  interactions 
involved. 

In  Ref.  7,  aerodynamic  constraints,  represented  by  induced  drag  constraints  are  added 
to  the  previous  set  of  constraints  explored.  New  approximations  for  induced  drag  constraints 
are  developed.  Again  the  composite  wing  of  an  RPV  is  used  for  numerical  experimentation 
with  the  new  capability.  Design  studies  with  design  variables  and  constraints  that  span  the 
disciplines  of  structures,  controls  and  aerodynamics  are  presented.  It  is  demonstrated  that 
the  synthesis  of  actively  controlled,  fiber  composite  wings  with  modeling  accuracy  that  is 
acceptable  for  preliminary  design  is  both  feasible  and  practical.  Results  of  the  new  design 
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studies  provide  interesting  and  important  new  insights  into  the  complex  nature  of 
multidisciplinary  interactions  present  in  wing  design. 

Reference  5  is  Dr.  Livne's  landmark  Ph.D.  dissertation  which  incorporates  both  the 
analysis  and  the  results  discussed  so  far;  in  a  single  document.  It  is  our  intent  to  publish  this 
document  as  a  Air  Force  Technical  Report  with  the  assistance  of  the  funding  agency. 

Recognizing  that  the  description  of  the  research  activity  provided  above  might  be 
fragmented,  we  find  it  useful  to  provide,  below,  a  concise  overview  of  the  principal  research 
accomplishments  described  in  Refs.  1-7.  This  research  has  produced  a  truly  integrated 
comprehensive  multidisciplinary  wing  synthesis  capability  which  is  the  most  advanced  of  its 
kind  available  to  date.  This  capability  which  is  depicted  schematically  in  Fig.  1,  integrates 
the  disciplines  of  structures,  aerodynamics,  and  controls.  There  is  no  comparable  capability 
available  in  the  academic,  industrial,  or  research  community. 

The  principal  accomplishments  are: 

1.  Fully  coupled  structural,  aerodynamic  and  control  system  analyses. 

2.  Extended  design  space,  where  one  can  change  simultaneously  structural,  aerod)  namic 
and  control  type  design  variables. 

3.  Inclusion  of  a  rich  mix  of  behavior  constraints  and  alternative  objective  functions. 

4.  Combination  of  analysis/sensitivity  capability  with  approximation  concepts  and 
optimization  algorithm. 

5.  Ability  to  deal  with  both  maneuver  and  handling  type  of  constraints. 
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6.  Demonstration  of  the  effectiveness  of  the  new  synthesis  procedure  by  applying  it  to 
optimization  of  conventional  wings,  swept  back  wings  and  forward  swept  wings  (RPV, 
F-16,  X-29). 

The  unique  features  of  this  synthesis  capability  are:  ■ 

1.  Balanced  high  quality,  yet  computationally  efficient  analysis  modeling. 

2.  Analytic  sensitivity  for  all  constraints  with  respect  to  all  design  variables. 

3.  It  contains  the  first  successful  approximations  of  aeroservoelastic  stability  and  gust 
response  constraints  in  the  design  space. 

4.  The  design  space  spans  three  disciplines  simultaneously. 

5.  The  formulation  includes  a  more  comprehensive  set  of  constraints,  design  variables 
and  alternative  objective  functions  than  structural  synthesis  programs  such  as  TSO, 
FASTOP  or  ASTROS. 

The  effectiveness  and  computational  efficiency  of  this  synthesis  capability  when 
applied  to  the  multidisciplinary  optimization  of  a  lightweight  fighter  is  depicted  in  Fig.  2. 
For  this  case  the  computing  times,  number  of  design  variables,  and  constraints  are  shown 
in  Fig.  2.  Such  a  problem  has  40  design  variables  and  over  3000  constraints.  The  combined 
computer  time  for  one  analysis,  and  sensitivity  analysis  is  about  6  min  of  CPU  time.  The 
number  of  analyses  for  a  converged  optimal  design  N*  is  between  10  and  20.  Thus  multidis- 
plinary  optimization  of  practical  airplane  fighter  wings  is  feasible  because  it  requires 
between  one  to  two  hours  of  CPU  time. 

The  synthesis  capability  which  has  been  described  above,  is  limited  to  the  subsonic 
flow  regime.  Another  important  contribution  made  in  the  course  of  this  research  contract 
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was  the  extension  of  the  aeroservoelastic  analysis  capability  to  the  transonic  flow  regime. 
Obviously  such  an  analysis  capability,  which  is  computationally  efficient,  is  a  prerequisite  for 
extending  this  synthesis  capability  to  the  transonic  flow  regime,  in  which  all  modern  fighter 
aircraft  must  operate. 

The  description  of  a  transonic  adaptive  aeroservoelastic  analysis  capability  based 
approximate  unsteady  time  domain  aerodynamics  is  presented  in  Refs.  8  and  9. 

References  8  and  9  describe  the  development  of  a  digital  adaptive  controller  capable 
of  suppressing  flutter  in  composi  te  wings  under  time  varying  flight  conditions  in  subsonic  and 
transonic  flow.  The  wing  structure  is  modeled  using  the  modeling  capability  such  as  that 
described  in  Ref.  1.  A  new  transonic  unsteady  aerodynamic  approximation  methodology  is 
developed  which  is  computationally  efficient  and  particularly  suitable  for  transonic 
aeroservoelastic  applications.  This  approximation  is  based  on  a  combination  of  unsteady 
subsonic  aerodynamics  with  a  transonic  correction  procedure.  The  transient  response  of  the 
aeroservoelastic  system  is  obtained  using  Roger’s  approximation  for  converting  frequency 
domain  aerodynamics  into  the  time  domain,  state  transition  matrices,  and  an  iterative  time 
marching  algorithm.  The  aeroservoelastic  system  in  the  time  domain  is  modeled  using  a 
deterministic  ARMA  model  together  with  a  parameter  estimator.  This  approach  enables 
one  to  compute  the  aeroelastic  flutter  boundaries  in  the  time  domain.  This  analysis 
capability  was  compared  with  experimental  data,  obtained  at  NASA  Langley,  and  good 
agreement  with  experimental  data  was  noted  for  the  low  transonic  Mach  number  range. 

The  linear  quadratic  controller  gain,  for  the  digital  flutter  suppression  system,  at  each 
time  step  is  obtained  using  an  iterative  Riccati  solver.  The  digital  adaptive  controller  is 
robust  with  respect  to  the  unknown  external  loads.  Flutter  and  divergence  instabilities  are 
suppressed  simultaneously  using  a  trailing-edge  control  surface  and  displacement  sensing. 
Acceleration  sensing  alone  is  inadequate  for  suppressing  static  instabilities  such  as 
divergence. 
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The  most  important  accomplishments  of  this  research  [8,9]  were: 

1.  This  study  represents  the  first  application  of  adaptive  optimal  control  methodology 
to  active  flutter  suppressions  in  transonic  flow. 

2.  Development  of  a  computationally  efficient  aeroservoelastic  analysis  capability  in 
transonic  flow  which  is  suitable  for  incorporation  in  a  general  integrated 
multidisciplinary  wing  design  synthesis  capability  such  as  that  described  in  Refs.  2 
through  7. 

During  the  course  of  this  research  contract  two  graduate  students  (Dr.  Eli  Livne  and 
Dr.  C.  Pak),  who  were  fully  supported  by  the  contract,  have  received  their  Ph.D.  degrees. 

It  is  evident  from  the  forgoing  that  this  research  contract  has  been  quite  fruitful.  In 
addition  to  the  substantial  productivity  represented  by  published  research  results,  this 
program  has  produced  truly  significant  and  original  contributions  to  the  state  of  the  art.  It 
has  also  provided  intellectual  stimulation  and  educational  enrichment  for  the  graduate 
students  affiliated  with  the  program.  It  is  our  earnest  hope  that  the  new  ideas  that  have 
grown  out  of  this  research  activity  will  significantly  influence  future  design  practice  in  the 
aerospace  industry. 
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RUN  TIME  FOR  LIGHT  WEIGHT  FIGHTER 
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APPENDIX  A 

Detailed  Summary  of  the  Research  Conducted 


Note:  This  Appendix  contains  Refs.  3, 4,  6,  and  8  cited  in  the  body  of  the  report 
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Towards  Integrated  Multidisciplinary  Synthesis  of  Actively 
Controlled  Fiber  Composite  Wings 


E.  Livne,*  L.  A.  Schmit.t  and  P.  P.  Friedmannt 
University  of  California,  Los  Angeles,  Los  Angeles,  California  90024 


The -synthesis,  of  nctlvdy  controlled  composite  wings  b  formulated  u  a  mulltdbdpllnxry  optimisation 
problem.  A  unique  Intefratlon  of  analysis  techniques  spanning  the  disciplines  of  structures,  aerodynamics,  and 
controls  is  described.  A  rich  variety  of  behavior  constraints  can  be  treated  Including  stress,  displacement,  control 
surface  travel  and  hinge  moment,  natural  frequency,  aeroservoelutic  stability,  gust  response,  and  handling 
quality  constraints,  as  well  as  performance  measures  la  terms  of  drag/lift  coefficients,  drag  polar  shape, 
required  load  factor  or  roll  rate,  and  wing  mass.  The  design  space  Indudes  a  simultaneous  treatment  of 
structural,  aerodynamic,  and  control  system  design  variables.  The  paper  sets  the  stage  for  multidisciplinary  wing 
optimization  by  describing  the  capabilities  and  discussing  the  accuracy  of  the  analysis  and  related  behavior 
sensitivity  analysis.  Applicability  of  approximation  concepts  to  the  mullidlKlplInary  optimization  problem  b 
examined  by  studying  typical  aeroservoelaslic  stability,  gust  response,  and  performance>related  constraints.  The 
computational  efficiency  of  the  combined  analysis  and  sensitivity  as  well,  as  the  quality  of  key  behavior 
constraint  approximations  Indicate  that  single-level  optimization  of  composite,  actively  controlled  practical 
wings  b  within  reach. 
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Nomenclature 

=  coefficients  of  the  numerator  polynomial 
in  a  sensor  transfer  function  (see  Fig.  3) 
=  Linear  time  invariant  (LT1)  state  space 
equation  matrices 
=  state  space  system  matrices  in 
standard  form 
=  aerodynamic  lag  terms 
=  aerodynamic  gust  lag  terms 
-control  law  transfer  function  numerator 
coefficients 

=  coefficients  of  the  denominator 
polynomial  in  a  sensor  transfer  function 
(see  Fig.  3) 

=  damping  matrix 

=  aerodynamically  modified  damping 
matrix  in  aeroservoelaslic  analysis 
=  transfer  function  denominator 
coefficients 

=  coefficients  of  the  numerator  polynomial 
in  an  actuator  transfer  function  (see 
Fig.  3) 

=  aeroelastic  system  matrices 
=  coefficients  of  the  denominator 
polynomial  in  an  actuator  transfer 
function  (see  Fig.  3) 

=  objective  function 

=  vector  of  behavior  constraint  functions 
=  aileron  gain  (see  Fig.  7) 

-stiffness  matrix 

=  aerodynamically  modified  stiffness 
matrix  in  maneuver  load  analysis 
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=  aerodynamically  modified  stiffness 
matrix  in  aeroservoelaslic  analysis 
=  mass  matrix 

=  aerodynamically  modified  mass  matrix  in 
aeroservoelastic  analysis 
=  mean  square  value  of  response  to  gust 

*  reference  mean  square  value  of  gust 
response 

=  number  of  actuator  states 
=  number  of  aerodynamic  added  states 
=  number  of  control  law  states 
=  number  of  control  surface  degrees  of 
freedom  (DOF) 

*  number  of  gust  filter  states 
=  number  of  structural  DOF 
=  number  of  sensor  states 

=  number  of  aeroservoelastic  states 
=  number  of  Roger  lag  terms 
=  number  of  Roger  lag  terms  added  for  gust 
=  any  design  variable 
=  load  vector 

=  minimum  state  or  Roger  approximation 
matrices 

=  generalized  displacements 
=  control  surface  deflection 
=  dynamic  pressure 
=  the  vector  of  Laplace  transformed 
generalized  displacements 
=  the  Laplace  transformed  gust  vector 
=  generalized  aerodynamic  force  matrix 
=  white  noise  intensity  matrix 
=  aerodynamic  added  states 
=  aerodynamic  gust  added  states 
=  Laplace  variable 
= reference  area 
=  skin  thickness  term 

=  vector  of  control  surface  excitations  (see 
Eq.  (39)) 

=  flight  speed 
=  white  noise  input 
=  vertical  gust  velocity 
=  the  Laplace  transformed  vertical  gust 
velocity 
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it/imm. 

■  aeroservoelastic  system  matrices 

w 

estate  vector 

X 

=  vector  of  design  variables 

M 

=  state  covariance  matrix 

y 

■  output 

t>, 

■  input  command  to  an  actuator 

=  structural  response  transformation 
matrices 

t 

■  equivalent  viscous  damping 

fo 

■  reference  equivalent  viscous  damping 

Subscripts 

0 

■  reference  value 

ACT 

■  actuator 

CO 

■  control  law 

c 

■  control  surface 

C 

■  gust 

s 

■  structural  DOF 

SE 

■  sensor 

Superscripts 

c 

■  gust 

L 

■  lower  bound 

V 

=  upper  bound 

aeroelastic  stability  design  problem  is  given  in  Ref.  76.  Para¬ 
metric  studies  reported  in  Refs.  7.’  and  78  further  highlight  the 
importance  of  multidisciplinary  design  considerations  for 
practical  fighter  wings.  An  integrated  approach  to  the  opti¬ 
mization  of  airplane  configuration  and  control  system  was 
presented  in  Refs.  79  and  80  using  simplified  mathematical 
models.  However,  the  application  of  modem  optimization 
techniques  to  wing  design  involving  a  diverse  mix  of  con¬ 
straints  based  on  analyses  from  several  disciplines  (structures, 
structural  dynamics,  aeroelasticity,  aerodynamics,  control, 
handling  qualities)  has  not  yet  been  treated  in  a  comprehensive 
and  realistic  manner. 

The  purpose  of  this  paper  is  *o  outline  a  unified  framework 
for  multidisciplinary  wing  synthesis.  It  contains  a  description 
of  a  set  of  techniques  for  analysis  and  behavior  sensitivity 
analysis  of  actively  controlled  fiber  composite  wings,  which 
will  facilitate  integrated  design  optimization.  Several  aspects 
of  structural,  aerodynamic,  and  control  system  modeling  are 
discussed.  Computational  efficiency  and  accuracy  of  approxi¬ 
mations  for  a  rich  variety  of  behavior  constraints  are  exam¬ 
ined.  The  paper  lays  the  foundation  for  the  application  of 
approximation  concepts  and  optimization  techniques  to  prac¬ 
tical  control  augmented  aeroelastic  wing  design. 


Introduction 


THE  introduction  of  active  control  technology1*4  and  com¬ 
posite  structural  tailoring1*"  to  airplane  wing  design  dur¬ 
ing  the  last  15  years  requires  a  re-examination  of  the  design 
practice  followed  in  the  past,  which  was  based  on  a  sequential, 
compartmented  approach. 

Using  the  sequential  approach,  undesirable  interactions  be¬ 
tween  disciplines,  which  were  not  taken  into  account  properly 
during  wing  development,  have  resulted  in  expensive,  some¬ 
times  lengthy,  modifications  and  fixes.12*16  At  the  same  time 
there  has  been  a  growing  recognition  of  the  potential  improve¬ 
ments  possible  when  an  integrated  multidisciplinary  design 
and  synthesis  approach  is  followed,  in  which  all  relevant  disci¬ 
plines  are  considered  simultaneously. ,7*19 
Three  decades  of  extensive  research  and  development  have 
made  optimization  techniques  widely  accepted  in  every  major 
discipline  needed  in  wing  design  synthesis.20  Structural  synthe¬ 
sis  has  matured  in  the  last  decade  and  has  already  been  used  to 
size  complex  structures  under  various  static  and  dynamic  be¬ 
havior  constraints  including  flutter  and  static  aeroelastic  con- 
straims.21*38  Research  carried  out  during  the  last  20  years  on 
control  system  design,  active  flutter  suppression,  and  gust 
alleviation  has  produced  several  alternative  implementations 
of  optimization  to  active  control  synthesis  practices.39-46  In  the 
aerodynamic  field,  optimization  has  been  used  to  synthesize 
wing  camber,  cross  section,  and  planform  to  achieve  desirable 
aerodynamic  characteristics  and  performance.47*31 
The  growing  confidence  in  modern  analysis  techniques  and 
disciplinary  optimization  has  led  to  a  departure  from  conven¬ 
tional  designs.  It  has  become  apparent  that  multidisciplinary 
interactions  must  be  taken  into  account  to  prevent  failure  and 
to  extract  maximum  benefits  from  the  design  freedom  offered 
by  a  truly  integrated  approach.  During  the  last  decade,  control 
augmented  structural  synthesis  has  emerged  as  an  important 
research  area.39*66  Initial  results,  for  wing  design,  have  been 
recently  reported  emphasizing  the  multidisciplinary  struc¬ 
tural/aerodynamic  synthesis  of  wings.67*71  An  initial  study  of 
the  aeroservoelastic  optimization  problem,  namely  the  simul¬ 
taneous  synthesis  of  wing  structures  and  their  active  control 
systems,  was  also  reported  recently.72  Methods  for  control  law 
and  control  system  performance  sensitivity  analysis  with  re¬ 
spect  to  structural  and  control  system  parameters  have  also 
been  reported  in  recent  years  as  a  precursor  to  the  application 
of  multilevel  decomposition  techniques  for  design  optimiza¬ 
tion.73*73  An  approach  to  the  integrated  handling  qualities/ 


The  Complexity  and  Multidisciplinary 
Nature  of  Wing  Design 

The  set  of  wing  design  descriptors,  whose  elements  consist 
of  preassigned  parameters  and  design  variables,22  is  shown  in 
Fig.  1.  Discussion  is  limited  to  wings  operating  in  the  sub- 
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sonic  to  low-supersonic  flight  speeds,  so  that  thermal  effects 
can  be  neglected.  A  given  set  of  design  descriptors  completely 
defines  a  particular  wing  design.  Which  of  the  descriptors  will 
be  preassigned  parameters  and  which  will  be  used  as  design 
variables  depends  on  the  level  of  application  for  optimization 
techniques  in  the  hierarchy  described  in  Ref.  22,  namely, 
whether  the  design  space  includes  sizing,  configuration  (ge¬ 
ometry),  or  topological  design  variables.  The  set  of  behavior 
functions,  from  which  constraints  aud  objectives  will  be  se¬ 
lected,  can  be  divided  into  two  categories  (Fig.  2).  Primary 
(system  level)  behavior  functions  are  those  performance  mea¬ 
sures  that  determine  the  overall  quality  and  competitiveness  of 
the  wing.  Secondary  (subsystem  level)  behavior  functions  are 
the  behavior  functions  that  must  be  taken  into  account  during 
the  design  to  guarantee  the  prevention  of  failure  in  all  possible 
failure  modes  and  to  introduce  known  constraints  on  subsys¬ 
tem  performance.  These  are  the  means  necessary  to  achieve 
the  overall  design  goals  and  should  ideally  be  “transparent" 
compared  with  real  design  objectives,  although  sometimes 
there  can  be  strong  correlation  between  a  secondary  behavior 
and  a  primary  behavior  function  (e.g.,  mass  and  airplane 
performance). 

The  importance  of  multidisciplinary  interactions  in  wing 
design  is  evident  from  Figs.  1  and  2.  Structural  topology, 
shape  and  sizing;  control  system  topology,  control  law  trans¬ 
fer  function  order,  and  gain  values;  as  well  as  aerodynamic 
configuration  layout,  jig  shape,  and  control  surface  deflec¬ 
tions  in  maneuvers  all  interact  to  achieve  desired  wing  perfor¬ 
mance  while  ensuring  structural  integrity,  aeroservoelastic  sta¬ 
bility,  ride  comfort,  and  good  handling  qualities.*2'93 

The  design  synthesis  problem  can  be  cast  in  a  mathematical 
programming  form: 

min, *,/■([*))  (1) 

s.t.  Igls(O) 

Uf‘|sWs|^l 

To  overcome  the  inherent  complexity  and  to  address  the 
computationally  intensive  nature  of  this  problem,  two  ap¬ 
proaches  have  been  suggested  in  the  literature.  The  first  ap¬ 
proach  is  based  on  the  application  of  multilevel  decomposition 
techniques  -ombined  with  existing  tools  for  detailed  analysis 
and  sensitivity  analysis  for  each  of  the  disciplines.94-93  The 
second  approach  seeks  to  gain  some  insight  into  the  nature  of 
the  problem  by  using  highly  simplified  mathematical  models 
or  simple  airplane  configurations  for  structural,  aerodynamic, 
and  control  system  analysis.72-79-*0 

Modeling  Considerations 

Structural  and  Aerodynamic  Modeling 

In  Ref.  20,  Ashley  points  out  that  a  considerable  part  of  the 
research  done  on  wing  optimization  has  been  based  on  models 
that  are  "a  long  way  from  the  complicated,  built  up  lifting 
surfaces  of  real  aircraft  with  their  multiple  design  criteria  and 
constraints.”  He  warns  that  “very  undesirable  consequences 
can  result  from  the  omission  or  careless  handling  of  con¬ 
straints,"  and  this  warning  is  particularly  relevant  for  mul¬ 
tidisciplinary  synthesis,  where  only  limited  experience  exists 
to  guide  the  designer,  and  intuition  may  sometimes  be 
misleading. 

The  prevalent  structural  beam/aerodynamic  strip  models 
used  for  basic  research  in  aeroelasticity  are  often  inadequate 
when  it  comes  to  synthesizing  real  wings.9*-9*  More  realistic 
models  are  needed.  At  the  same  time,  detailed  finite  element 
models  and  Computational  Fluid  Dynamics  (CFD)  aerody¬ 
namic  techniques  are  still  computationally  too  expensive  to 
use  within  the  inner  loop  of  a  multidisciplinary  synthesis  ap¬ 
proach.  Thus,  it  is  necessary  to  bridge  the  gap  between  the 
highly  idealized  and  the  very  detailed  modeling  alternatives  by 


introducing  balanced  analysis  and  design  optimization  models 
that  capture  essential  behavior  characteristics,  without  making 
the  integrated  multidisciplinary  design  optimization  task  in¬ 
tractable. 

The  integrated  optimum  design  capability  outlined  here  is 
based  on  modeling  and  analysis  techniques  for  the  required 
disciplines,  which  are  consistent  with  each  other  in  terms  of 
accuracy  and  efficiency  and,  thus,  lead  to  a^balanced  treat¬ 
ment  of  practical  wings.  In  the  structures  area,  a  rather  gen¬ 
eral  equivalent  plate  analysis,9*  which  builds  on  the  basic  ideas 
underlying  the  Aerodynamic  Tailoring  and  Structural  Opti¬ 
mization-  (TSO)  computer  code7-2*-29;  and  incorporates  addi¬ 
tional  recent  developments  due  to  Giles, W-IC0  is  used.  The 
equivalent  plate  approach  for  structural  modeling  of  low  as¬ 
pect  ratio  wings  has  been  known  for  many  years.  It  was  Giles, 
however,  who  showed  that,  using  present  day  computers,  a 
single  high-order  power  series  can  be  used  for  approximating 
displacements  over  wing  planforms  made  of  several  trape¬ 
zoidal  segments  to  obtain  accurate  stress  as  well  as  displace¬ 
ment  information.  Stresses  in  spar  and  rib  caps  can  be  calcu¬ 
lated  in  addition  to  composite  skin  stresses.  Configurations 
made  of  several  plate  segments  attached  to  each  other  via 
springs  accounting  for  attachment  stiffness  and  actuator  stiff¬ 
ness  can  be  analyzed  to  simulate  wing/control  surface  config¬ 
urations.  The  simplicity  of  manipulating  simple  power  series 
leads  to  analytic  rather  than  numerical  integ.ation  for  the 
mass  and  stiffness  expressions.  With  the  careful  organization 
of  computer  storage  space  and  ordering  of  calculations,  major 
savings  in  computation  times  and  core  storage  requirements 
can  be  achieved. 

In  the  work  described  here,  the  equivalent  plate  structural 
analysis  documented  in  Ref.  98  is  integrated  with  the  Piece- 
wise  Continuous  Kernel  Function  Method  (PCKFM)  devel¬ 
oped  by  Nissim  and  Lottati  for  lifting  surface  unsteady  aero¬ 
dynamics.101'104  Lifting  surface  unsteady  aerodynamics103-106 
has  served  as  the  basic  aerodynamic  modeling  tool  for  the 
flutter  analysis  of  airplanes  since  the  1960s.  The  PCKFM 
combines  the  power  of  the  doublet  lattice  method  in  dealing 
with  pressure  singularities  with  the  accuracy  and  speed  of  the 
kernel  function  method.  Extensive  numerical  experimentation 
has  demonstrated101  that  PCKFM  is  accurate  and  converges 
rapidly.  For  configurations  involving  control  surfaces,  it  can 
take  narrow  gaps  into  account,  is  faster  than  lattice  methods, 
and  is  more  accurate  in  the  calculation  of  control  surface 
hinge  moments.  Thus,  it  is  particularly  suitable  for  calculating 
the  generalized  unsteady  air  loads  (on  lifting  surfaces  made  up 
of  wing  and  control  surface  elements)  that  are  needed  for 
active  flutter  suppression  and  gust  alleviation  studies. 

The  combination  of  modern  equivalent  plate  structural 
modeling  and  PCKFM  lifting  surface  aerodynamics  is  thought 
to  be  adequate  for  the  preliminary  design  of  airplane  wings 
and  for  the  exploratory  venture  into  multidisciplinary  practi¬ 
cal  wing  synthesis.  In  addition  to  a  reliable  prediction  of 
flutter  results  and  static  aeroelastic  effects,  useful  hinge  mo¬ 
ment106  and  induced  drag  predictions16’-10*  can  be  expected  for 
subsonic  and  supersonic  small  angle-of-attack  flight.  The 
analysis  is  adequate  for  addressing  flight  stability  and  control 
problems  of  the  elastic  airplane.109  Its  aerodynamic  predic¬ 
tions  might  be  improved  by  using  correction  factor  techniques 
if  any  measured  data  are  available. 

Ftalte-Dimeftsioiia]  State  Space  Modeling  of  Uaslrady  Aerodynamics 

It  is  suggested  in  the  literature  that  aeroservoelastic  stability 
analysis  can  be  successfully  based  on  the  p-k  method  using 
generalized  aerodynamic  force  matrices  computed  for  simple 
harmonic  motion.37-110  However,  when  the  optimization  of  the 
design  for  aeroservoelastic  stability  is  addressed  and  modem 
control  techniques  are  to  be  implemented,  it  is  necessary  to 
cast  the  aeroelastic  equations  of  motion  in  Linear  Time  Invan- 
ant  (LTI)  state  space  form.  It  then  follows  that  some  approx¬ 
imation  of  the  unsteady  aerodynamic  loads  in  terms  of  ratio¬ 
nal  functions  of  the  Laplace  variable  is  needed. 
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Thf  method  of  Roger111  has  been  widely  used  for  finite-di¬ 
mensional  unsteady  aerodynamic  loads  representation  during 
the  past  decade.  A  series  of  aerodynamic  stiffness,  damping, 
apparent  mass,  and  several  lag  terms  is  used  to  approximate 
elements  of  the  generalized  aerodynamic  loads  and  gust  force 
matrices  over  a  range  of  reduced  frequencies.  A  least-squares 
fitting  procedure  is  carried  out  for  each  element  of  the  ma¬ 
trices  separately. 

The  minimum  state  method,  developed  by  Karpelm  and 
recently  studied  in  Ref.  113,  is  found  to  be  attractive  because 
it  has  the  potential  for  generating  accurate  approximations  to 
unsteady  generalized  aerodynamic  forces,  while  adding  only  a 
small  number  of  states  to  the  mathematical  model  of  the 
aeroservoelastic  system.  It  is  based  on  an  iterative  fitting  pro¬ 
cess  in  which  all  terms  of  the  generalized  aerodynamic  load 
and  gust  force  matrices  are  considered  simultaneously.  In 
comparison  with  other  finite  state  modeling  techniques,  the 
number  of  states  needed  in  the  minimum  state  method  appears 
to  be  smaller  for  the  same  overall  accuracy  of  approxima¬ 
tion.113  This  leads  to  a  state  space  model  of  lower  order,  thus 
reducing  core  requirements  and  computation  time. 

Considering  the  accumulated  experience  and  fast  generation 
of  Roger  approximations  (resulting,  however,  in  higher-order 
mathematical  models  of  the  aeroservoelastic  system)  vs  the 
smaller-order  mathematical  models  possible  with  the  mini¬ 
mum  state  approach  (with  a  relative  lack  of  experience  and 
time-consuming  approximant  generation  as  potential  handi¬ 
caps),  it  was  decided  to  include  both  methods  in  the  present 
capability  as  available  alternatives. 


Control  System  Modeling 

The  integrated  aeroservoelastic  system  is  modeled  as  an  LTI 
system.  Since  the  number  of  sensors  and  control  surfaces  is 
small  in  real  airplanes,  the  complex,  high-order  laws  generated 
by  some  multivariable  control  system  design  techniques  are 
avoided  at  this  stage.  A  schematic  block  diagram  of  the  ac¬ 
tively  controlled  aeroservoelastic  system  is  shown  in  Fig.  3. 
Airplane  motions  (acceleration  and  angular  rates)  are  mea¬ 
sured  by  a  set  of  sensors  placed  on  the  structure.  The  resulting 
signals  are  used  as  inputs  to  the  control  law  block  which 
commands  control  surface  actuators.  The  control  surface  mo¬ 
tions  guarantee  stability  and  desirable  dynamic  response  of  the 
complete  system. 

The  control  system  is  completely  described  by  locations  of 
sensors  and  control  surfaces  and  by  the  transfer  functions  of 
the  sensors,  control  laws,  and  actuators.  Gain  scheduling  can 
be  adopted  by  assigning  different  control  laws  to  different 
flight  conditions. 


Optimization  Considerations 
Design  Variables 

Shape  design  variables  have  already  received  considerable 
attention  in  wing  optimization  studies.67'7l•7,•,0',5  However,  in 
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Fig  3  Elastic  airplane  control  system  block  diagram. 


addition  to  a  balanced  approach  to  analysis  in  terms  of  the 
analysis  techniques  selected,  we  seek  to  keep  a  balance  in  level 
of  optimization  by  focusing  at  present  on  sizing-type  design 
variables  in  all  disciplines  considered  (Fig.  1  j. 

Structural  Design  Variables 

Figure  4  shows  an  airplane  modeled  as  an  assembly  of 
flexible  lifting  surfaces.  Each  lifting  surface  is  modeled  as  an 
equivalent  plate  whose  stiffness  is  controlled  by  contributions 
from  thin  cover  skins  (fiber  composite  lamintes)  and  the  inter¬ 
nal  structure  (spar  and  rib  caps).  Wing  sections  are  connected 
to  each  other  via  stiff  springs  (representing  hinge  stiffness  at 
attach  points)  and  flexible  springs  (representing  the  stiffness 
of  actuators  and  their  backup  structure).  Each  wing  section 
can  include  several  trapezoidal  parts.  Concentrated  masses  are 
used  to  model  nonstructural  items  and  balance  masses. 

The  vertical  displacement  wof  each  wing  section  is  approx¬ 
imated  by  a  Ritz  polynomial  series  of  the  form 

w  {x,y,t)~  £  <7,  (2) 

!•] 

where  x  and  y  are  chordwise  and  spanwise  coordinates  respec¬ 
tively.  The  exponents  m,  and  n,  define  the  specific  polynomial 
series  used.  It  can  be  a  complete  polynomial  in  x  and  y  or  a 
product  of  polynomials  in  x  and  y  (see  Ref.  98). 

The  depth  of  a  wing  section  is  given  by  a  polynomial 

A'a 

h  (x,y)=  £  H,  x’yu  (3) 

i-  ■ 

where  the  H,  are  preassigned  parameters. 

Thickness  distribution  of  a  typical  skin  layer  is  represented 
by 

Kx,y)=t  T,xly'-  (4) 

i- 1 

Rib  and  spar  cap  areas  are  allowed  to  vary  linearly  along  their 
length 


A[t])  =  A0  +  A,  i)  (5) 

Wing  stiffness  and  mass  matrix  elements  are  linear  combina¬ 
tions  of  certain  area  integrals  over  the  planform,  line  integrals 
over  spar/rib  length,  and  polynomial  terms  evaluated  at 
points  where  concentrated  masses  are  located  or  springs  are 
attached. M  The  present  equivalent  plate  modeling  capability’* 
makes  it  possible  to  efficiently  analyze  combined  wing  box/ 
control  surface  configurations.  A  wing  assembly  and  a  canard 
or  horizontal  tail  may  be  attached  to  a  fuselage  (modeled  as  a 
flexible  beam  or  a  flexible  plate)  to  simulate  complete  airplane 
configurations.  The  level  of  modeling  detail  can  be  selected 
independently  for  each  section.  Therefore,  the  degree  of  detail 
used  to  model  control  surfaces  for  analysis  and  synthesis  is  not 
limited,  as  is  the  case  of  the  TSO  code. 

At  the  present  stage  of  research,  structural  topology,  shape, 
and  materia]  properties  are  preassigned;  however,  skin-layer 
fiber  orientations  are  available  as  design  variables.  For  sktn- 


Fig.  4  Airplane  *s  an  nsjeasbly  of  tqalvnieai  plate. 
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layer  thicknesses  [Eq.  (4))  the  coefficients  of  the  thickness 
power  series  serve  as  design  variables.  This  guarantees  smooth 
thickness  variation  for  each  layer.  For  spar  and  rib  cap  areas 
[Eq.  (5)1,  two  coefficients  are  used  as  de  gn  variables  for  each 
spar  or  rib.  Concentrated  masses  at  preassigned  locations  and 
spring  constants  for  linear  and  rotational  springs  can  also  be 
treated  as  design  variables. 


Aerodynamic  and  Control  System  Design  Variables 

Wing  cross  section,  aerodynamic  planfurm,  and  topology 
are=preassigned  here.  Performance  and  loads  in  quasistatic 
maneuvers  can  be  influenced  by  designing  the  jig  shape  (initial 
camber)  of  the  wing  and  by  proper  deflection  ofleading-edge 
and  trailing-edge  control  surfaces.  The  initial  camber  of  the 
wing  is  given  by  a  series 

n . 

»°(x,y,t)  =  £  <7,°(r)  xmyni  (6) 

i-i 

where  the  powers  m,  and  n,  are  identical  to  those  in  Eq.  (2), 
and  any  subset  of  the  coefficients  <7®  can  serve  as  design 
variables.  The  deflections  of  control  surfaces  for  each  distinct 
maneuver  point  are  also  available  as  separate  design  variables. 

The  control  system  design  variables  at  the  lowest  level  in  the 
hierarchy  (analogous  to  sizing)  are  the  coefficients  of  numera¬ 
tor  and  denominatoriof  control  law  transfer  functions.  Con¬ 
trol  surface  locations,  sensor  locations,  topology  of  the  con¬ 
trol  system,  and  order  of  numerator  and  denominator 
polynomials  in  the  transfer  functions  are  preassigned.  It  is  also 
assumed  that  sensor  and  actuator  transfer  functions  are  preas¬ 
signed,  although  the  formulation  is  sufficiently  general  so  as 
to  allow  the  treatment  of  their  numerator  and  denominator 
coefficients  as  additional  design  variables. 

The  set  of  design  variables  treated  spans  three  disciplines, 
namely  structures,  aerodynamics,  and  control.  The  design 
space  is  thus  opened  up  to  include  sizing  level  design  variables 
from  all  three  disciplines  simultaneously. 

Behavior  Functions 

In  order  to  provide  for  a  rich  variety  of  constraints  and 
alternative  objective  functions,  the  following  analysis  capabil¬ 
ities  are  included. 

1)  Static  ‘‘maneuver  load"  analysis  (static  aeroelastic  de¬ 
flection  and  stress  calculations  for  the  elastic  airplane  in  ma¬ 
neuver).  maneuvers  include  symmetric  pull-ups  (defined  by 
Mach  number,  altitude,  and  load  factor)  or  steady  rolling 
maneuvers  (defined  by  Mach  number,  altitude,  and  roll  rate). 
In  addition  to  elastic  deflections  and  stresses,  the  control 
surface  deflections  and  hinge  moments  needed  for  the  maneu¬ 
vers  are  calculated. 

2)  Static  ‘‘given  loads"  analysis  (static  deflection  and  stress 
calculations  for  the  cantilevered  wing  under  a  set  of  prescribed 
loads),  the  loads  are  assumed  independent  of  the  structural 
design  and  do  not  change  in  the  course  of  wing  synthesis.  This 
option  is  important  for  cases  where  linear  aerodynamic  theory 
is  inadequate,  forcing  the  use  of  experimental  data. 

3)  Natural  frequency  and  mode  shape  analysis,  natural  fre¬ 
quencies  and  mode  shapes  are  obtained  for  different  sets  of 
boundary  conditions  (this  facilitates  generation  of  separate 
symmcric  or  antisymmetric  modes). 

4)  Aero.-rvoelastic  stability  analysis,  poles  of  the  control- 
augmented  an,.lane  are  calculated  for  different  level  flight 
conditions  (defined  by  specifying  Mach  number  and  altitude). 

5)  Gust  response  analysis,  root-mean-square  (rms)  values  of 
control  surface  rotations  and  rates  and  rms  values  of  selected 
sensor  measurements  due  to  continuous  atmospheric  turbu¬ 
lence  are  calculated  for  different  flight  conditions. 

6)  Drag  analysis,  induced  drag  is  calculated  for  the  elastic 
lift  distribution  during  maneuvers.  Drag  values  assuming  ei¬ 
ther  full  leading-edge  suction  (fully  attached  flow)  or  no  lead¬ 


ing-edge  suction  (separated  flow  at  the  leading  edge)  or  a 
combination  of  these29,35  are  available. 

With  the  integrated  analysis  capability  that  has  been  deve.. 
oped,  the  following  behavior  functions  can  be  evaluated,  elas 
tic  displacements;  elastic  twist,  spar/rib  cap  stresses,  sk.n 
combined  stress  failure  criteria"4;  natural  frequencies,  real 
and  imaginary  parts  of  aeroservoelastic  poles;  rms  values  of 
random  control. surface  rotations  and  rates  due  to  gust;  rms 
values  of  sensor  measurements  in  gust;  total  mass;  lift  and 
drag  coefficients;  control  surface  rotations;  and  hinge  mo¬ 
ments  in  maneuvers.  It  should  also  be  noted  that  roll  rate  or 
load  factor  at  a  particular  altitude  and  Mach  number  can  be 
treated  as  design  variables;  therefore,  they  can  be  maximized 
or  constrained. 

Control  surface  effectiveness  is  not  addressed  directly  at  this 
stage.  The  synthesis  emphasizes  sustaining  a  desired  roll  rate 
or  load  factor  while  keeping  hinge  moments,  control  surface 
deflections,  and  stresses  within  allowable  bounds. 

Aeroservoelastic  stability  is  guaranteed  by  providing  ade¬ 
quate  damping  at  each  flutter  critical  aeroservoelastic  pole 
throughout  the-flight  envelope.115  Handling  qualities  can  be 
preliminarily  addressed  via  inequality  constraints  on  the  aero¬ 
servoelastic  pole  locations  (e.g.,  short-period  root  placement) 
and  pilot-seat  acceleration  due  to  atmospheric  turbulence. 4419 
The  control  surface  deflection  needed  for  trim  and  overall 
performance  in  a  given  maneuver  and  its  rms  activity  due  to 
gusts  can  be  combined  to  ensure  that  no  saturation  occurs.9,1 9: 

Any  of  the  behavior  functions  or  their  combinations  can 
serve  as  objective  functions.  Possible  alternatives  are  mass, 
drag  (to  be  minimized),  steady  roil  rate,  lift-to-drag  ratio  (to 
be  maximized),  or  a  combination  of  these. 

The  present  analysis  capability  offers  a  rich  variety  of  be¬ 
havior  functions  for  wing  design  synthesis  Thus,  the  interac¬ 
tion  among  structure,  control,  aerodynamics,  handling  quali¬ 
ties,  and  airplane  performance  can  be  taken  into  account  in  an 
integrated  manner. 

Approtch  to  Integrated  Optimization 

Once  the  preassigned  parameters,  design  variables,  failure 
modes,  load  conditions,  and  objective  function  are  selected, 
the  integrated  optimization  problem  can  be  cast  as  a  nonlinear 
programming  problem  having  the  form  of  Eq  (1) 

The  nonlinear  programming  approach  combined  with  ap¬ 
proximation  concepts  (NLP/AC  approach)  has  proven  to  be 
an  effective  method  for  solving  structura’  synthesis  prob¬ 
lems,21 22  and  here  it  will  be  adapted  to  the  m.’tidisciplinary 
design  optimization  task.  In  this  method  relatively  few  de¬ 
tailed  analyses  are  carried  out  during  optimization.  Each  anal¬ 
ysis  and  the  associated  behavior  sensitivity  analysis  serve  as  a 
basis  for  constructing  approximations  to  the  objective  and 
constraint  functions  in  terms  of  the  design  variables.  Thus,  a 
series  of  explicit  approximate  optimization  problems  is  solved 
converging  to  an  optimal  design. 

The  main  advantage  of  the  NLP  formulation  is  its  general¬ 
ity.  No  a  priori  assumptions  have  to  be  made  about  the  set  of 
active  constraints  at  the  optimum  Given  an  initial  design,  a 
local  optimum  is  sought  using  mathematical  programming 
techniques.  Thus,  it  is  especially  suitable  for  multidisciplinary 
optimization,  where  the  problem  is  large  and  complicated  and 
past  experience  does  not  provide  much  intuitive  guidance 
However,  for  the  NLP  approach  to  be  practical,  it  is  crucial  to 
avoid  too  many  detailed  analyses.  Success  in  this  regard  de¬ 
pends  on  efficient  analysis/sensitivity  calculations  and  oh 
making  the  explicit  approximations  of  objective  and  con¬ 
straint  functions  robust  yet  simple  enough  for  efficient  solu¬ 
tion. 

The  use  of  analytic  behavior  sensitivities  and  the  construe 
tion  of  robust  approximations  for  behavior  functions  in  term, 
of  the  design  variables  are  at  the  heart  of  the  NLP  AC  ap¬ 
proach.  Luring  the  past  two  decades,  approximation  tech 
mques  for  static  deflection,  stress,  and  natural  frequency  eon 
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straints  have  been  studied  extensively  and  are  now  well 
established.  Several  methods  for  divergence  and  flutter  con¬ 
straint  treatment  have  beendeveloped,  but  no  experience  has 
been  reported  with  the  aeroservoelastic  poles  of  an  actively 
controlled  airplane,  the  rms  of  the  response  to  random  gusts, 
or  hinge  moments  of  a  flexible  control  surface  on  an  elastic 
wing.  It  should  be  pointed  out  that  in  maneuver  load  analysis 
the  loads  acting  on  the  airplane  are  functions  of  the  design 
variables  via  static  aeroelasticity  as  well  as  through  inertial 
effects.  The  right  side  |P)  of  the  matrix  equation 

(Alt?  )  =  IP)  (7) 

(where  k  is  a  stiffness  matrix  modified  by  aerodynamic  terms 
corresponding  to  structural  and  control  surface  motions)  thus 
depends  on  the  structural  design  variables,  whereas  in  the 
classical  given  loads  analysis,  the  right  side  l  P )  is  fixed.  Thus, 
the  success  of  approximations  using  reciprocal  variables  (Ref. 
116)  for  stress  and  static  deflection  in  static  analysis  is  not 
guaranteed  in  maneuver  load  analysis,  and  alternative  approx¬ 
imations  have  to  be  carefully  examined. 

Some  aspects  of  maneuver  load  and  drag  analysis  have 
already  been  discussed  in  the  literature  within  the  framework 
of  integrated  w  ing  optimization.29  33'6’  64  10  In  the  following  we 
will,  therefore,  focus  on  the  aeroservoelastic  and  gust  response 
analysis  and  behavior  sensitivity  calculations. 


Aeroservoelastic  Analysis 

As  pointed  out  in  Ref.  117,  all  time-dependent  phenomena 
of  the  elastic  airplane  are  governed  by  a  universal  set  of 
equations  of  motion,  wherein  only  the  right  side  (representing 
the  input)  varies  in  proceeding  from  one  phenomenon  to  the 
next.  Indeed  a  measure  of  multidisciplinary  integration  of 
analysis  techniques,  concepts,  and  terminology  is  needed  even 
before  multidisciplinary  optimization  is  addressed. 

Several  steps  in  this  direction  have  been  taken  in  the 
past  n*-i2i  However,  almost  20  years  after  the  publication  of 
Ref.  118,  still  no  particular  approach  to  the  dynamics  of  the 
deformable  airplane  is  universally  accepted.  The  analysis  here 
for  time-dependent  problems  is  based  on  the  widely  used  set  of 
linearized  equations  for  small  perturbations  of  elastic  airplane 
motion  with  respect  to  constant  speed,  level  flight. ,22'127  Al¬ 
though  perturbations  in  the  longitudinal  direction  are  not 
taken  into  account  at  this  stage,  the  analysis  is  adequate  for 
addressing  basic  stability  and  control  as  well  as  aeroservoelas¬ 
tic  problems  of  highly  augmented  airplanes  in  the  context  of 
preliminary  design.*9 

AeroeUstlc  Model 

Assuming  small  perturbations  from  a  steady  level  flight,  the 
Laplace  transformed  equations  of  motion  for  the  elastic  air¬ 
plane  are 

|  [A/)s2+[C]s +  (*]  )(?(*)) 

-<7z>SlQ(5)](9(s)!  =?z>S(e0(J)!  (8) 

The  vector  of  displacements  can  be  expressed  in  terms  of  n , 
structural  response  DOF  and  n(  control  surface  deflections: 


and  the  equations  of  motion  corresponding  to  the  structural 
DOF  are  partitioned  accordingly  to  yield 

I  [Mu]s^[Cd]s  +  [A*])  [q,  |  +  f  [A/„]s2  +  (CW Js 
+  \ku] )  \q, )  -<7dS[Qd(*)]  1  I?j  I 

-qoS\QAs)))\q<)=*j-S  ICc,(s)I  wc(*>  (10) 


The  finite-dimensional  approximation  to  the  Laplace  trans¬ 
formed  generalized  unsteady  loads  using  the  RogeT  ap- 
proach111,113  is  of  the  form 


(Q(s)]=s2Ip,)+s(p2]+[p3]+,E -4r;lp'->)  (>» 

In  the  minimum  state  method,"2  "3  the  functional  dependence 
of  the  generalized  aerodynamic  force  matrix  on  the  Laplace 
variable  is  approximated  by  a  rational  expression  of  the  form 

[Q(s)l  =  |P,1*}  +  IP&  +  iPy)  +  \P*)[sV)  ~  IPs)]  *  'IPJ*  (12) 

The  finite-dimensional  Roger  approximation  of  the  Laplace 
transformed  generalized  unsteady  gust  loads  used  is  of  the 
form 


\Qc(s))=s[Pf]  +  \P?\  +  I P?.  i) 


(13) 


The  minimum  state  approximation  used  for  the  gust  load 
vector  is 

[Qg(s)\  =  [Pf)s  + 1 P?)  +  \P?\  [*!/)  -  |FSC)]  *  '[P6c)s  (14) 

(Note  that  the  notation  IP)  is  used  to  denote  matrices  associ¬ 
ated  with  unsteady  aerodynamics  finite-dimensional  state 
space  approximations.  However,  the  matrices  [P,]  and  |P,C) 
(/&  4)  have  different  meanings  depending  on  whether  the 
Roger  or  minimum  state  approximation  is  used.) 

Partitioning  the  aerodynamic  matrices  associated  with 
Structural  DOF  according  to  Eq.  (9)  leads  to  a  Roger  approx¬ 
imation  of  the  form 

[Q”(s).e,r(5)l  =  s2[P^PH+ ♦  [Pf.Pf] 


/•A'/  £ 

+  (?I  (s  +  bi) 


[P*  j.  Piiy  1 


(15) 


and  to  a  minimum  state  partitioned  approximation  of  the 
form 

tQ°(s)l  =  tPf]s2  +  [Pf)5  +  [Pf]  f  [Pills/  -  P5)  -  'IPils  (16) 


\Qu(s))  =  [P,*>:+  [Pf]s  +  IPfl  +  (Pills/ -P«)-'  (Pels  (17) 

In  the  above  expressions,  IPf),[Pf),|Pfl  ate  aerodynamic 
apparent  mass,  damping,  and  stiffness  matrices,  respectively. 
Their  dimensions  are  (n,xn,).  lP/t],lPjr).|Pf)  are  aerody¬ 
namic  apparent  mass,  damping,  and  stiffness  matrices,  respec¬ 
tively,  associated  with  the  coupling  terms  between  control 
surface  and  structural  motion.  Their  dimensions  are  (n,  xne). 

Aerodynamic  states  are  now  introduced  as  follows.  For  the 
Roger  approximation 


(18) 

(19) 

Aerodynamic  states  for  the  minimum  state  method  are  in¬ 
troduced  as 

itfWI  =  [sJ/J-IPslJ-'lP.’ Pfr  [*£j] 

(20) 

|/f(s))  =  [s(/]-  [Pfl]  - '  (Pf)  H-c  (S ) 

CD 
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Thus,  in  the  Roger  approximation,  the  state  space  equations 
for  the  added  aerodynamic  states  associated  with  structural 
deflection  are 

+  (22) 

and  in  the  minimum  state  approximation 

J I  r*  (s) )  =  (Pj)  f^(j))+  (Pfls  I  q,(s) )  +  [P^s  ( qc(s) )  (23) 

Two  vertical  wind  gust  state  space  models  are  available  in 
the  literature  for  random  gust  response  calculations.  These  are 
the  Dryden5-4  model  and  any  rational  approximation  to  the 
von-Karman,ajN  model.  In  both  cases  a  gust  filter,  repre¬ 
sented  by  a  strictly  proper  transfer  function,  is  used  to  trans¬ 
form  a  Gaussian  zero  mean  white  noise  input  w  into  the 
vertical  gust  speed  wc  with  a  given  power  spectral  density  and 
rms.  A  state  space  description  of  the  gust  input  is  thus  of  the 
form 


4(xc  1  =  [-4c]lxc]  +  (Be  )w 

(24) 

*c(4)  =  [Cc](Jfc] 

(25) 

*  »c(s)=lCc)lAcUxc )  +  [Cc)[Bc  1  w 

(26) 

If  aerodynamicaily  augmented  mass,  damping,  and  stiffness 
matrices  are  now  introduced  in  the  form 

(27) 

Ca  =  C„-qDSP2a 

(28) 

X^Kn-qoSP? 

(29) 

M^^Mx-qnSP? 

(30) 

C«  =  Ctc-qDSP? 

(31) 

X^Kx-qoSP? 

(32) 

(Kk  and  CK  are  zero) 

and  five  subvectors  of  the  state  vector  (x)  are  defined  as 

(*.)  =  (?,! 

(33) 

(x2  ]  =s|?,| 

(34) 

(xj  ]  =  [xq  1 

(35) 

* 

II 

T 

ik 

(36) 

(*5)  =  |rc! 

(37) 

then  after  substitution  of  Eqs.  (11-32)  into  Eq.  (10),  the 
structural  and  unsteady  aerodynamic  part  of  the  integrated 
equations  of  motion  can  be  written  in  matrix  form  as 


*f£)i*)=ini*)+[G]{u)+i/ntv 

where  |u  ]  is  a  vector  of  control  surface  excitations 


\u\={  sq< 

(aj 


(38) 


(39) 


while  J£),(F],(G],  and  [H]  are  matrices  whose  elements  de¬ 
pend  on  structural,  aerodynamic,  and  gust  filter  terms.  Equa¬ 
tions  (38)  are  the  Laplace  transformed  equations  for  the  struc¬ 
tural  and  aerodynamic  states.  To  complete  the  state  space 


formulation  of  this  problem,  the  feedback  expression  relating 
the  control  surface  motion  |&(s))  to  the  structural  motions 
(?i(i)li  Uftts)),  UJft(*))  is  used  to  close  the  control  loops. 

Control  System 

The  actual  displacements,  velocities,  and  accelerations  at  a 
set  of  points  on  the  structure  are  given  by 

bsTRfr)!  =  (  [*o]  +  l*ll*  +  [*2k2|f<7.(*))  (40) 

The  matrices  *0.  +i.  *2  are  determined  by  the  location  of 
measurement  points  on  the  structure.  Actual  structural  re¬ 
sponses  are  measured  by  a  set  of  sensors  whose  output  signals 
serve  as  input  to  controlflaws.  These  control  laws  generate 
input  commands  (6)  to  the  actuators.  The  n,  commands,  5,, 
are  synthesized  to  serve  as  input  to  nc  actuators. 

A  typical  control  law  transfer  function  is  given  by 

&i  _(2ffs"  +  b„-,s"-|  + . +  b0 

^seT  + . +  d„ 

where  >Sei  is  the  input  signal  to  the  control  element.  The  state 
equations  for  a  single  control  law  are  given  by 


s  (Xcol  = 


0  0 
1  0 
0  1 


0  0 


■do 

-d, 

-di 


-4,-i 


l*co) 


(42) 
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6/  =  (0  0  0  •  •  •  —  1 )  |*co)  +  4,  -^se/ 


(43) 


where  d,  and  b,  are  the  transfer  function  denominator  and 
numerator  coefficients.  When  all  control  laws  arc  assembled, 
the  multi-input  multi-output  control  block  controlling  several 
actuators  is  ’■"presented  by 

■r!*co)  =  (Aco)l*co)  +  [ScoII^se)  f44) 


|6)  =  (Q:o!l*co)  +  (AtoII-Vse)  (45)  | 

Sensor  and  actuator  transfer  functions  are  assumed  to  have  1 

denominators  of  higher  order  than  their  numerators.  When 
several  sensors  are  present,  the  state  space  model  relating  the 
sensor  states  |xSe).  the  sensor  outputs  |.Vse).  and  the  actual 
structural  response  |.ySTR )  is 

*  ( *SE  )  =  MseI  1*SE  1  +  (5se1  I ^STR  I  (46) 

(>'se)  =  ICse](xSe)  (47) 


The  state  space  model  relating  the  sensor  states  lxSE  | .  the 
sensor  outputs  |y$E),  and  the  actual  response  |>strI  >5 
Eqs.  (33),  (34),  (40),  and  (47)]. 
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s  Use!  =  MseJUse )  +  [5se)([^o) (■*! ) 

+  [*i]l*z)  +  [+z*)(jrj))  (48) 

The  state  space  model  of  the  actuators  is  given  by 

s  UactI  =  MactII^act)  +I5act)U!  (49) 

iQc )  =  (CactIUact)  (50) 

where  ( qc )  is  the  control  surface  deflection  (assumi  i  , 
versible  surfaces).  The  matrices  [ACo],  [ficol.  [Ax„,  1 
(5se],  Mact).  and  [Fact!  depend  on  the  des:gn  var1'1 '  “  ?  .  j 
ciated  with  the  control  system  through  reiatienr 
Eqs.  (42)  and  (43).  Using  Eqs.  (44-50),  the  co;:-v-  sur'.-..-- 
excitation  vector  («)  (see  Eq.  (39)]  can  now  be  erF.  ,-sr.vd  .n 
terms  of  structural,  aerodynamic,  and  control  system  sta-'.s  as 
follows: 

(?<■  1  =  [CactI  !*act1  «■  S ) 

5 \Qc  1  =  (CactIMact) Uact!  +  [CactH-SactIICco] Ucc  * 

+  (CAcr)[5ACT][/>co][CsE]  (*se  )  (52) 

1  =  [CactHBactII^coHCseIIBse! 

X (l*o)  1  Qs  I  +  (4>|]s  (9,)  +  [4>j]s J ( q , ) ) 

+  ICact!(-4act]:  Uact  1  +  [Cact]  ([/4  actH^act]  l-Dcol !  Cse] 
+  ISactIICcoM^coHCse]  +  [Sact][Dco][Cse](/1se]) 
x  Use!  +  ICAcr!(MACr][5AcT](Cc0] 

+  [£act][Qto]I'4co])Uco!  (53) 

Complete  System 

The  state  vector  of  the  whole  sysit  t  is  partitioned  into  eight 
subvectors.  The  following  sub\ect>  -  \re  added  to  the  five 
which  have  been  already  defined  in  Eqs.  (33-37). 


l*6l  =  l*AC7)(»ACrXl) 

(54) 

(*l)  =  I*seI(”seX  1) 

(55) 

l*s)  =  l*col(ncoXl) 

(56) 

Assembly  of  the  state  space  models  of  the  structure,  sensors, 
actuators,  control  block,  and  approximate  unsteady  load  and 
gust  aerodynamics  leads  to  the  closed-loop  state  space  equa¬ 
tions  of  the  complete  system: 

*[(/)(*(*)}  =  [»'][*(*)}  +  (»')>*’ 

(57) 

The  order  of  the  system  is  given  by 

zisys  =  2  n,  +  nACT  +  nSE  +  nc0  +  «c  +  «aer 

(58) 

When  N,  lag  terms  are  used,  then  in  the  case  of  Roger  approx¬ 
imation 

«aer  -N,nt  +  Nf 

and  in  the  case  of  minimum  state  approximation 
^aer  -Nf+Np 

The  considerable  saving  in  terms  of  added  aerodynamic  states 
with  the  minimum  state  approach  is  now  evident,  since  with 
the  Roger  approximation  .he  number  of  added  states  is  a 


product  of  the  number  of  lag  terms  useJ  to  fit  the  data  [Eq, 
(11)]  times  the  number  of  structural  DOF  used. 

Dependence  of  System  Murker  on  tie  Design  Varieties 

Examine=the  [£],  [F],  [G],  and:[/f)  matrices  of  Eq.  (38). 
When  augmented  by  the  control  system  state  space  equations 
for  the  control  system  states  [Eqs.  (44),  (48),  and  (49)]  and  if 
full-order  structural  matrices  are  used  (no  modal  reduction), 
then  each  of  these  matrices  is  a  linear  combination  of  struc¬ 
tural,  aerodynamic,  and  products  of  control  system  terms. 
The  transformation  matrix  relating  system  states  |x  |  to  con¬ 
trol  excitations  [«)  [Eqs.  (51-53))  depends  only  on  control 
system  parameters  (no  modal  reduction).  Substitution  of  Eqs. 
(51-53)  into  Eq.  (38)  shows  that  the  elemen.s  of  the  [U],  (H. 
and  [  W ]  matrices  of  Eq.  (57)  are  of  the  form 

(SI),,  +  (A  1)„  +  (Cl)„  +  (S2)„(C2)„  +  (A  2)„(C3)„  (59) 

The  structural  mass  and  stiffness  terms  SI,  S2  depend  only- 
on  structural  design  variables,  whereas  the  control  system 
matrices  Cl,  C2,  C3  depend  only  on  control  system  design 
variables  via  the  state  space  models  of  actuators,  sensors,  and 
control  laws.  Since  wing  cross  section  and  planform  shape  are 
preassigned  at  present,  when  full-order  stiffness  and  mass 
matrices  are  used,  the  aerodynamic  terms  associated  with  gen¬ 
eralized  unsteady  loads  and  gust  forces  do  not  change  with  a 
change  in  design  variables.  When  moda  reduction  is  used  to 
reduce  the  order  of  the  system  matrices,  then  a  fixed  modes 
approach  is  adopted  here1*0  resulting  in  fixed  aerodynamic 
terms  for  the  modally  reduced  models  also.  Modal  reduction  is 
further  discussed  in  the  next  section. 

Full-Order  and  Reduced-Order  Models 

With  the  equivalent  plate  approach,  the  structural  model 
may  include  a  relatively  small  number  of  DOF  when  compared 
to  conventional  finite  element  models.  These  are  generalized 
displacements  associated  with  the  Ritz  polynomials.  The  gen¬ 
eralized  aerodynamic  matrices  are  calculated  for  the  same  set 
of  Ritz  polynomials  used  in  Eq.  (2).  Thus,  the  aeroservoelastic 
stability  analysis  can  be  done  with  a  full  order  model,  includ¬ 
ing  the  full-order  mass,  stiffness,  and  aerodynamic  matrices  or 
by  order  reduction  based  on  a  small  number  of  normal  modes. 

If  a  Roger  or  a  minimum  state  approximation  can  be  found 
that  will  accurately  fit  the  full-order  aerodynamic  matrices 
with  a  small  number  of  lag  terms,  then  these  approximation 
matrices  can  be  used  in  the  aeroservoelastic  stability  analysis 
even  with  modal  reduction.  They  just  have  to  be  prcmultiplied 
end  postmultiplied  by  the  generalized  mode  shapes  in  order  to 
have  the  approximation  to  the  modally  reduced  aerodynamic 
matrices.  Although  the  number  of  structural  DOF  might  dif¬ 
fer  significantly  between  full-order  and  modal  analysis,  the 
number  of  aerodynamic  states  is  the  same  when  both  ap¬ 
proaches  are  based  on  the  same  minimum  state  approximation 
to  the  full-order  aerodynamic  matrices.  If  the  Roger  approxi¬ 
mation  is  used,  the  number  of  added  aerodynamic  state,  in¬ 
creases  with  any  increase  in  the  number  of  modes  used  for 
reduction.  This  produces  very  large  models  when  many  modes 
or  a  full-order  analysis  are  used.  Another  possibility  is  first  to 
modally  reduce  the  frequency-dependent  unsteady  aerody¬ 
namic  matrices  using  modes  that  are  periodically  updated 
after  a  given  number  of  gualysis/ssniitivitv  optimization  cy¬ 
cles.  Then  it  might  be  possible  to  fit  these  reduced  matrices 
using  fewer  lag  terms  than  the  number  needed  to  fit  the 
full-order  matrices.  If  a  smaller  number  of  lag  terms  can  thus 
be  used,  this  will  reduce  the  dimensions  of  the  U,  V  matrices 
in  Eq.  (57),  compared  to  their  dimensionality  in  a  modal 
approach,  which  uses  p/egenerated  full-order  aerodynamic 
matrices.  However,  rather  than  preparing  the  full-order  aero- 
Gynamic  approximantes  only  once  before  any  synthesis  stans, 
in  the  latter  case  it  would  be  necessary  to  generate  aerody¬ 
namic  approximants  each  time  the  set  of  modes  used  as  a  basis 


for  modal  reduction  is  changed  during  the  synthesis.  All  of  the 
modal  reduction/unsteady  aerodynamic  approximation  alter¬ 
natives  described  so  far  have  been  included  in  the  present 
capability.  However,  their  comparative  assessment  is  beyond 
the  scope  of  this  paper. 

When  full-order  aerodynamic  matrices  are  used,  they  are 
generated  once  for  the  Ritz  functions  employed  in  the  struc¬ 
tural  analysis,  and  they  are  invariant  with  respect  to  changes  in 
either  structural  sizing  or  control- system  design  variables. 
When  modal  reduction  is  used,  the  mode  shapes  are  periodi¬ 
cally  updated  after  a  given  number  of  analysis  steps  and  so  are 
the  aerodynamic  matrices  and  their  finite-dimensional  rational 
approximation:.  Nevertheless,  in  eigenvalue  sensitivity  calcu¬ 
lations,  the  derivatives  of  all  reduced-order  matrices  are  deter¬ 
mined  using  a  fixed-mode  approach.150  This  might  require  the 
use  of  more  modes  in  order  to  obtain  good  sensitivity  informa¬ 
tion.  In  summary,  the  derivatives  ofithe  aerodynamic  matrices 
with  respect  to  structural  or  control  system  design  variables 
are  zero  for  the  full-order  case,  and  they  are  assumed  zero  for 
the  reduced-order  case. 

The  derivative  aU/dp,  dV/dp  can  be  calculated  analyti¬ 
cally  for  either  structural  or  control  systenrudesiiiiv  variables 
(p)  by  differentiting  Eq.  (59).  The  gust  load  vector  {W\ 
depends  only  on  gust  filter  and  aerodynamic  design  variables, 
and  therefore  its  partial  derivative  with  respect  to  p  vanishes. 

Stability  is  examined  by  computing  the  eigenvalues  of  the 
generalized  eigenvalue  problem: 

Mt/(P)JI*1-l»'(P>]|*)  (60) 

Sensitivity  of  eigenvalues  with  respect  to  structural  and  con¬ 
trol  system  design  variables  is  given  by 


where  d  and  «  are  left  and  right  eigenvectors,  respectively. 

Gust  Response  Annlysb  and  Sensitivity 
In  addition  to  several  publications  addressing  it  in  the  con¬ 
text  of  active  control  technology, 3-{-39  airplane  response  to 
random  atmospheric  turbulence  has  already  been  discussed  in 
the  context  of  structural  optimization. ,3I>132  Here,  attention  is 
focused  on  the  rms  values  of  control  surface  rotations  \qe ), 
rates  (qc).  and  sensor  measurements  (ysEl- 


The  state  space  equations  [Eq.  (57)]  are  transformed  into 
standard  fonn 

S  (x(s))*[X)(x(s))  +  (5)w(j)  (62) 

Since  only  (qc).  r(<?t  1.  and  lysEl  are  considered,  [Eqs.  (47), 
(51),  and  (52)],  it  follows  that  each  output  considered  yt  is 
given  by 

>t  =  |C*!rU)  (63) 

where  (C*  j  is  either  constant  or  a  function  of  control  system 
design  variables. 

The  state  covariance  matrix  is  a  solution  of  a  Lyapunov’s 
matrix  equation135  in  the  form 

i/fim+m(^F*-[^)[g.][^)r  («) 


where  [Q„]  is  the  intensity  matrix  of  th;  Gaussian  white  noise 
w.  Sensitivity  of  the  covariance  matrix  [Sf]  with  respect  to  a 
design  variable  p  is  calculated  by  differentiating  Eq.  (64): 


and  finally  the  derivative  of  the  (rms)2  of  yk  is  calculated  by 
differentiation  of  the  covariance  expressions  based  on  Eqs. 
(47),  (51),  and  (52).  It  should  be  noted  that 


i/n=if/]-'[K] 

(66) 
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(67) 

and 


from  which  it  follows  that 


Fig.  5  Lightweight  fighter  structure]  model. 


Fig.  6  Lightweight  fighter  ocrodyMmlc  model. 
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Fig.  7  Lightweight  fighter  control  system  roll  cb«nnel. 
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Fig.  S  Approximations  or  missile  pitch  root  dimping  (variable  wing 
thickness). 


and 


The  Hessenberg-Schur  method  of  Ref.  134  is  used  to  solve 
both  Eqs.  (64)  and_(65).  Note  that  after  the  analysis  is  carried 
out,  the  matrices  [A ),  [A]r  are  already  reduced  to  Hessenberg 
and  Schur  forms.  Therefore,  the  sensitivity  calculations  of  F.q. 
(65)  are  equivalent  to  adding  right  sides  to  Eq.  (64). 

Numerical  Examples 

Actively  Controlled  Lightweight  Fighter  Teat  Model 
Structural  and  aerodynamic  models  of  a  lightweight  fighter 
airplane  are  shown  in  Figs.  5  and  6.  The  airplane  is  similar  to 
the  YF-16,  and  the  construction  of  its  mathematical  model 
was  guided  by  Refs.  29,  135,  and  136.  It  is  different.than  an 
actual  YF-16.  Vertical  tail  and  ventral  fins  are  not  included. 
The  fuselage  and  horizontal  tail  are  assumed  rigid.  There  is  no 
leading-edge  flap  The  wing  has  a  biconvex  (about  49o  thick- 
ness-to-chord  ratio)  airfoil.  It  is  made  of  aluminum  skins  and 
an  array  of  spars  and  ribs.  A  flexible  aileron  and  a  rigid  tip 
launcher/missile  assembly  are  attached  to  the  wing  using 
springs  The  configuration  analyzed  weighes  20,000  lb  and  is 
statically  unstable. 

The  roll  channel  of  the  flight  control  system  is  shown  in  Fig. 
7.  It  is  based  on  the  YF-16  roll  channel  (Refs.  12,  135,  and 
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Fig.  9  Approximations  of  missile  pitch  root  damping  (variable  roll 
loop  gain). 


136).  The  lightweight  fighter  model  used  here  is  intended  to 
illustrate  the  well-known  YF-16  aeroservoelastic  roll  instabil¬ 
ity  (Refs.  12,  13,  and  135).  Lack  of  sufficient  data  precluded 
a  more  refined  simulation  of  the  YF-16.  However,  the  model 
used  here  is  representative  of  a  typical  fighter  airplane  and  is 
quite  complex  in  its  detail. 

Two  mechanisms  of  instability ,  similar  to  those  encountered 
in  the  YF-16,  exist  here.  As  the  aileron  gain  Fp  (Fig.  7)  is 
increased,  a  6.5  Hz  instability  appears  associated  with  a  mis¬ 
sile  pitch  mode.  With  a  further  increase  of  Fp,  a  second 
instability  appears  at  3.5  Hz  associated  with  the  rigid-body  roll 
mode.  Mach  0.9  aerodynamics  for  antisymmetric  motion  at 
20,000  ft  is  used  in  the  stability  calculations.  Since  the  original 
airplane  is  unstable  at  this  point,  it  is  artifitially  stabilized  here 
by  reducing  the  gain  Fp  and  adding  inherent  damping.  Stabil¬ 
ity  of  the  model  is  necessary'  for  studying  gust  response  ap¬ 
proximations. 

Results 

The  multidisciplinary  wing  analysis/synthesis  capability  de¬ 
scribed  here  has  been  extensively  tested.  Structural,  aerody¬ 
namic,  aeroelastic,  and  aeroservoelastic  results  have  been 
compared  with  analysis/test  data  available  from  other 
sources,  and  overall  good  correlation  was  found.  The  analytic 
behavior  sensitivities  were  verified  by  finite-difference  sensi¬ 
tivity  checks.  Here  we  emphasize  some  of  the  basic  issues 
associated  with  the  feasibility  of  applying  the  NLP/ AC  ap¬ 
proach  to  truly  integrated  muludisciplinary  wing  synthesis. 
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These  are  the  computational  speed  of  analysis/sensitivity 
calculations  and  the  robustness  of  behavior  function  approxi¬ 
mations. 

Parametric  studies  of  the  effects  of  design  variable  varia¬ 
tions  on  aeroservoelastic  poles  are  presented  in  Figs.  8  and  9. 
Figure  8  shows  the  variation  of  the  damping  ratio  f  of  the 
missile  pitch  mode  with  a  typical  structural  thickness  design 
variable.  In  Fig.  9  the  damping  varies  with  a  typical  control 
system  design  variable,  the  gain  (Fp)  (Fig.  7).  In  both  figures 
the  design  variable  is  varied  over  a  wide  range  of  values.  In 
practical  optimization,  move  limits  will  be  imposed  on  design 
variables  to  ensure  accuracy  of  behavior  function  approxima¬ 
tions.  When  move  limits  of  30%  are  used,  both  direct  and 
reciprocal  Taylor  series  representations' 16  yield  reasonable  ap¬ 
proximations.  Thus,  hybrid  approximations  seem  adequate 
for  explicit  representations  of  aeroservoelastic  pole  con¬ 
straints. 

Figure  10  shows  the  variation  of  the  aileron  mean  square 
deflection  due  to  atmospheric  turbulence  when  the  gain  Fp  is 
varied.  As  the  missile  pitch  mode  is  destabilized  (Fig.  9),  the 
gust  response  increases  sharply.  Away  from  the  stability 
boundary,  the  mean  square  (ms)  aileron  deflection  is  well-ap¬ 
proximated  by  either  a  linear  or  reciprocal  approximation  with 
30%  move  limits.  Tighter  move  limits  might  be  needed  near 
the  stability  boundary  as  damping  ratios  approach  zero.  It  is 
expected  that  any  optimization  procedure  will  respond  to  a 
decrease  in  damping  and  an  increase  in  gust  response  by 
driving  the  design  away  from  damping  and  gust  response 
constraint  violations. 

Constraints  associated  with  maneuver  loads  are  evaluated 
next  to  determine  the  quality  of  approximations  (see  Eq.  (7)). 
The  thickness  of  the  skin  on  the  wing  inboard  box  is  varied 
Over  a  wide  range  of  values.  The  variation  of  the  aileron  hinge 
moment  needed  to  sustain  a  desired  steady-state  roll  at  sea 
level,  M  =  0.9,  is  shown  in  Fig.  11.  A  stress  constraint  for  a 
point  on  the  skin  in  a  9  g  symmetric  pull-up  in  terms  of  a 
quadratic  stress  failure  criterion"4  is  illustrated  in  Fig.  12.  The 


turbulence  (variable  roll  loop  gain). 


Fig,  II  Approximations  of  bilge  moment  in  steady  roll  (variable 
wing  ibkkness). 


aileron  hinge  moment  increases  as  aileron  effectiveness  is  kx, 
due  to  a  decrease  in  wing  skin  thickness  in  Fig  1 1 .  Aga  r.  with 
move  limits  of  30%,  both  linear  and  reciprocal  approxima. 
tions  work  quite  well.  Based  on  the  examples  presented  here, 
it  is  expected  that  robust  approximations  can  be  constructed 
for  constraints  associated  with  a)  deflection,  stress,  cor/ro! 
surface  trim  angle,  and  hinge  moment  in  given  pull  up  and  roll 
maneuvers  and  b)  aeroservoelastic  stability  and  gust  respc  se 
limitations  in  given  level  flight  conditions  Similar  considc 
ations  can  be  expected  to  apply  to  deflection  and  stress  co' 
straints  in  “given-load"  load  conditions  as  well  as  natura’ 
frequency  constraints  based  on  the  extensive  study  of  these 
constraints  within  the  framework  of  structural  synthesis 
Typical  CPU  times  for  execution  of  an  analysis  /sernsiti v  it  \ 
cycle  on  the  lightweight  fighter  are  given  in  Table  1  In  this 
example,  three  maneuver  load  cases  and  one  given-load  case, 
symmetric  and  antisymmetric  modes  are  calculated,  and 
aeroservoelastic  stability  and  gust  response  are  included  for 
one  level  flight  condition.  The  wing  is  covered  by  a  grid  of  81 
points  for  deflection  and  stress  calculations  A  totai  of  3222 
constraints  and  all  of  their  sensitivities  with  respect  to  40 
structural,  aerodynamic,  and  control  system  design  variables 
are  calculated  in  about  7  min.  The  constraints  include  a  com 
prehensive  mix  of  gauge,  slope,  displacement,  stress,  natural 
frequency,  aeroservoelastic  pole,  gust  response,  drag,  hinge 
moment,  and  control  surface  travel  constraints  As  Table  I 
shows,  the  major  part  of  the  computation  time  in  one  analy 
sis/sensitivity  cycle  is  spent  on  the  behavior  sensitivity  calcula 
tions.  Constrain!  deletion*2  strategies  will  reduce  this  time 
considerably  by  retaining  only  a  small  subset  of  critical  and 
potentially  critical  constraints  as  drivers  in  each  approximate 
problem  formulation.  Only  sensitivities  of  the  retained  con 
straints  are  needed,  and  CPU  time  for  one  detailed  analysis  ’ 
sensitivity/approximate  problem  generation  stage  will  be  re 
duced  considerably.  Thus,  if  about  10-15  detailed  analyses  are 


Fig.  12  Approximations  of  skin  quadratic  (tress  failure  criterion  in  a 
9  |  symmetric  pull-up  maneuver  (variable  wing  thickness). 


Table  1  Typical  computation  times,  seconds 


Generate  M,  K,  and  |P]  5.4 

Given-load  solution  0.09 

Maneuver  load  solutions  0.35 

Drag  calculations  0.0' 

Natural  modes  4.8 

Generate  A,  B  matrices  1 .53 

Aeroservoelastic  subility  analysis  15.55 

Gust  response  analysis  11.44 

Deflection,  stresses  calculations  0.85 

Total  analysis  40,12 

Stess,  deflection  sensitivities  295,15 

Natural  frequency  sensitivities  1.4“ 

Aeroservoelastic  pole  sensitivities  + 
gust  response  sensitivities  53.15 

Total  sensitivity  349.~~ 


.  niA'-rvw  i 


needed  for  optimization  based  on  approximation  concepts 
(Ref.  22),  it  can  be  anticipated  that  between  40-60  CPU  min 
will  be  needed  on  the  UCLA  IBM  3090  Model  200  for  inte¬ 
grated  multidisciplinary  optimization  of  practical  wings. 


Concluding  Remarks 

A  general  framework  for  wing  optimization  has  been  devel¬ 
oped,  highlighting  the  multidisciplinary  nature  of  the  prob¬ 
lem.  A  balanced  multidisciplinary  wing  analysis  and  behavior 
sensitivity  analysis  capability  has  been  described.  Emphasis 
was  placed  on  various  aspects  of  the  aeroservoelastic  problem 
formulation  as  well  as  integration  and  testing  of  the  aeroelas- 
tic  elements  of  the  new  method.  Promising  results  in  terms  of 
approximation  accuracy  and  computation  times  indicate  that 
the  integrated  multidisciplinary  optimization  of  practical  ac¬ 
tively  controlled,  fiber  composite  wings  is  within  reach. 
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Abstract 

Analysis  and  synthesis  techniques  used  in  a  newly 
developed  multidisciplinary  control  augmented  fiber 
composite  wing  optimization  capability  arc  reviewed. 
Structural,  aerodynamic  and  control  system  math¬ 
ematical  models  that  arc  suitable  for  the  preliminary 
design  of  real  airplanes  arc  used  in  an  integrated 
manner  to  synthesize  improved  designs  of  wings  and 
their  active  control  systems.  Optimization  techniques 
developed  for  structural  synthesis  arc  adapted  to  the 
integrated  multidi>ciplinary  wing  synthesis  problem, 
in  which  constraints  from  several  disciplines  arc 
taken  into  account  simultaneously  and  the  design 
space  is  opened  up  to  include  structural,  control 
system  and  aerodynamic  design  variables.  The  effec¬ 
tiveness  and  efficiency  of  the  new  capability  are 
studied  using  a  mathematical  model  of  a  remotely 
piloted  vehicle  (RPV). 

Introduction 

Multidisciplinary  interactions  have  always  been  at 
the  heart  of  airplane  wing  design.  The  introduction 
of  active  control  technology1*3  and  composite  struc¬ 
tural  tailoring1**  during  the  last  fifteen  years  have 
made  these  interactions  more  complex  and  more 
important-  Recent  experience  has  shown  that  not 
accounting  properly  for  multidisciplinary  interactions 
during  the  design  process  can  lead  to  dangerous 
consequences7'1.  At  the  same  time  the  benefits  of 
multidisciplinary  integration  have  become  widely 
recognized  motivating  extensive  research  and  influ¬ 
encing  design*' ". 
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Still  the  design  practice  of  the  past,  which  was 
based  on  a  sequential,  compartmcntcd  approach,  is 
followed  today.  True,  advanced  anahsis  and  testing 
tools  have  been  developed  to  address  multidiscipli¬ 
nary  interaction.  Active  flutter  suppression,  maneu¬ 
ver  load  control,  gust  alleviation  and  ride  smoothing 
techniques  arc  utilized,  to  mention  a  few  examples. 
Optimization  techniques  arc  used  for  control  system 
synthesis,  aerodynamic  design  and  aeroclastic  tailor¬ 
ing.  Structural  wing  synthesis  subject  to  structural, 
aeroclastic  and  aerodynamic  performance  constraints 
has  successfully  followed  in  the  footsteps  of  struc¬ 
tural  synthesis*1*13  using  a  variety  of  computer  codes 
and  tcchniqucsM*u. 

However,  the  application  of  optimization  tech¬ 
niques  is  still  done  one  discipline  at  a  time.  True 
integration  in  design,  namely  the  synthesis  of  wines 
subject  to  a  set  of  multidisciplinary  constraints 
addressing  design  variables  from  several  disciplines 
simultaneously,  has  not  been  carried  out. 

Especially  complex  and  difficult  is  the  integrated 
synthesis  of  wing  structure  and  its  active  control 
system.  In  the'  few  eases  where  integrated  synthesis 
was  studied'*'3*  the  mathematical  models  used  were 
so  simplified  that  they  did  not  provide  an  assessment 
of  the  techniques  needed  to  optimize  realistic  wings, 
and  the  accumulation  of  practical  design  experience 
was  not  possible.  Thus,  the  application  of  modem 
optimization  techniques  to  wing  design  involving  a 
diverse  mix  of  constraints  and  design  variables  based 
on  analyses  from  several  disciplines  (structures,  struc¬ 
tural  dynamics,  aeroclasticity,  control,  handling  qual¬ 
ities)  has  not  yet  been  treated  in  a  comprehensive 
and  realistic  manner. 

Reference  22  describes  3  synthesis  capability  for 
actively  controlled  fiber  composite  lifting  surfaces.  It 
discusses  structural,  aerodynamic  and  control  system 
modeling  techniques  adopted  and  modified  to  chal¬ 
lenge  the  formidable  computational  task  of  carrying 
out  analysis  and  behavior  sensitivity  calculations  for 
this  problem  with  computer  resources  that  will  make 
optimization  practical  and  with  accuracy  that  is  suffi¬ 
cient  for  preliminary  design.  It  is  based  on  the 
generalization  of  the  nonlinear  programming 
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approximation  concepts  (NLP/AC)12-'3  approach 
from  structural  synthesis  to  multidisciplinary  opti¬ 
mization. 

The  present  paper  reports  first  results  of  wing 
control  augmented  structural  synthesis  achieved 
using  the  new  capability.  The  applicability  of 
approximation  concepts  to  the  control  augmented 
structural:  synthesis  of  wings  is  studied.  Examination 
of  optimization  convergence  as  influenced  by  includ¬ 
ing  approximations  of  new  constraints  (especially 
aeroscrvoclastic  stability  and  gust  response 
constraints)  guides  identification  of  effective  move 
limits,  convergence  criteria  and  approximation  types 
to  be  used.  Examples  of  integrated  optimization  of 
realistic  v.ings  and  their  active  control  systems  with 
structural  and  control  system  design  variables  subject 
to  gage,  stress,  aeroscrvoclastic  stability  and  gust 
response  constraints  offer  an  improved  understanding 
of  this  complex  synthesis  problem. 

Review  of  Mathematical 
Modeling  Techniques 

Structural  Modeling 

The  integrated  optimum  design  capability  is 
based  on  approximate  analysis  techniques  for  the 
required  disciplines,  which  are  consistent  with  each 
other  in  terms  of  accuracy  and  efficiency  and  thus 
lead  to  a  balanced  treatment  of  practical  wings.  In 
the  structures  area,  a  rather  general  equivalent  plate 
analysis23,  vvnich  builds  on  the  basic  ideas  underlying 
the  TSO  computer  code14  and  incorporates  addi¬ 
tional  recent  developments  proposed  by  Giles24,  is 
used.  High  order  simple  power  series  are  used  for 
approximating  displacements  over  wing  planforms 
made  of  several  trapezoidal  segments  to  obtain  accu¬ 
rate  stress  as  well  as  displacement  information. 
Stresses  in  spar  and  rib  caps  can  be  calculated  in 
addition  to  composite  skin  stresses.  Configurations 
made  of  several  plate  segments  attached  to  each 
other  via  springs  (accounting  for  attachment  stiffness 
and  actuator  stiffness)  can  be  analyzed  to  simulate 
wing '‘control  surface  configurations. 

Aerodynamics 

The  equivalent  plate  structural  analysis  docu¬ 
mented  in  Ref.  23  is  integrated  with  the  PCKFM 
(Piecewise  Continuous  Kernel  Function  Method) 
developed  by  Nissim  and  Lottati  for  lifting  surface 
unsteady  aerodynamics23  2*.  This  method  is  partic¬ 
ularly  suitable  for  calculating  the  generalized 
unsteady  air  loads  (on  lifting  surfaces  made  up  of 
wing  and  control  surface  elements)  that  are  needed 
for  active  flutter  suppression  and  gust  alleviation 
studies. 


When  the  optimization  of  the  design  for.acroser- 
voelastic  stability  is  addressed  and  modem  control 
techniques  are  to  be  implemented,  it  is  necessary  to 
cast  the  acroelastic  equations  of  motion  in  I.TI 
(Linear  Time  Invariant)  state  space  form.  It  then 
follows  that  some  approximation  of  the  unsteady 
aerodynamic  loads  in  terms  of  rational  functions  of 
the  l  .aplace  variable  is  needed. 

The  method  of  Roger”  has  been  widely  used  for 
finite  dimensional  unsteady  aerodynamic  loads  repre¬ 
sentation  for  quite  a  while.  The  Minimum  State 
Method,  developed  by  Karpel30  and  recently  studied 
in  Ref,  31,  is  found  to  be  attractive  because  it  has 
the  potential  for  generating  accurate  approximations 
to  unsteady  generalized  aerodynamic  forces,  while 
adding  only  a  small  number  of  augmented  states  to 
the  mathematical  model  of  the  aerosenoelastic 
system.  Both  methods  are  available  in  the  present 
capability. 

Control  System 

The  integrated  aeroscrvoclastic  system  is  modeled 
as  a  Linear  Time  Invariant  (I.TI)  system.  Since  the 
number  of  sensors  and  control  surfaces  is  small  in 
real  airplanes,  the  complex,  high  order  control  laws 
generated  by  multivariable  control  system  design 
techniques  are  avoided  at  this  stage,  A  block 
diagram  of  the  actively  controlled  aeroservoclastic 
system  is  shown  in  Figure  1.  Airplane  motions 
(acceleration  and  angular  rates)  are  measured  by  a 
set  of  sensors  placed  on  the  structure.  The  resulting 
signals,  ySE,  are  used  as  inputs  to  the  control  law 
block  which  commands  control  surface  actuators 
'Hie  control  surface  motions,  qc,  guarantee  stability 
and  desirable  dynamic  response  of  the  complete 
system. 

The  control  system  is  completely  described  by 
the  location  of  sensors  and  control  surfaces  and  by 
the  transfer  functions  of  the  sensors, control  laws  and 
actuators.  Gain  scheduling  can  be  adopted  by  assign¬ 
in')'  different  control  laws  to  different  flight  condi¬ 
tions. 

Modeling  Capabilities 

The  combination  of  modem  equivalent  plate 
structural  modeling  and  I’CRTM  lifting  surface  aero¬ 
dynamics  is  assumed  to  be  adequate  for  this  explora¬ 
tory  venture  into  multidisciplinary  practical  wing 
synthesis.  In  addition  to  a  reliable  prediction  of 
deformations,  stresses,  flutter  results  and  static  acroe¬ 
lastic  effects,  quite  good  hinge  moments32  and 
induced  drag33-34  can  be  expected  for  subsonic  and 
supersonic  small  angle  of  attack  flight  The  analysis 
is  adequate  for  addressing  flight  stability  and  control 
problems  of  the  elastic  airplane33.  Its  aerodynamic 
predictions  might  be  improved  by  using  correction 
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factor  techniques  if  any  measured  data  are  available. 
As  to  control  system  modeling,  the  techniques  used 
here  make  it  possible  to  properly  model  real  flight 
control  or  active  flutter  suppression  terns,  Fffects 
of  real  actuators  and  sensors  on  r.’rol  system 
performance  arc  automatically  takci.  into  account 
through  their  given  transfer  functions. 

Design  Variables 

Ref.  22  sets  forth  a  framework  for  multidiscipli¬ 
nary  wing  synthesis  and  describes  a  hierarchy  of 
design  variables  consisting  of  sizing  type  design  vari¬ 
ables  at  the  lowest  level,  followed  by  shape  and  then 
toplogical  type  design'  variables,  This  classification 
applies  to  design  variables  spanning  struc¬ 
tures, control  and  aerodynamics.  Here,  in  addition  to 
a  balanced  approach  to  behavior  prediction  in  terms 
of  the  analysis  techniques  selected,  a  balance  is  main¬ 
tained  in  level  of  optimization  by  focusing  on  sizing 
type  design  variables  for  all  disciplines  considered. 

Structural  Design  Variables 

figure  2  shows  an  airplane  modeled  as  an  assem¬ 
bly  of  flexible  lifting  surfaces.  Fach  lifting  surface  is 
modeled  as  an  equivalent  plate  whose  stiffness  is 
controlled  by  contributions  from  thin  cover  skins 
(fiber  composite  laminates)  and  the  internal  structure 
(spar  and  rib  caps).  [Mate  sections  are  connected  to 
each  other  via  stiff  springs  (representing  hinge  stiff¬ 
ness  at  attach  points)  and  flexible  springs  (represent¬ 
ing  the  stiffness  of  actuators  and  their  backup 
structure).  Fach  wing  section  can  be  made  of  several 
trapezoidal  parts.  Concentrated  masses  are  used  to 
model  nonstructural  items  and  balance  masses. 

The  vertical  displacement,  w,  of  each  wing 
section  is  approximated  by  a  Ritz  polynomial  series 
of  the  form 

,v. 

=£<?,<()  *m'/'  (1) 

1 

"here  x  and  y  are  chordwisc  and  spanwisc  coordi¬ 
nates  respectively,  m,  and  n,  are  powers  reflecting  the 
t\pe  of  polynomial  scries  used.  It  can  be  a  complete 
polynomial  in  x  and  y  or  a  product  of  polynomials 
in  x  and  y  (Ref.  23).  The  depth  of  a  wing  section  is 
gixen  by  a  polynomial 
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where  the  //,  are  preassigned  parameters.  Thickness 
Ji'tnhution  of  a  typical  skin  layer  is  represented  by 
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Rib  and-spar  cap  areas  are  allowed  to  vary  linearly 
over  their  length  t] 

A(m)=Ao  +  A,>7  (4) 

The  present  equivalent  plate  modeling  capability21 
makes  it  possible  to  efficiently  analyse  combined 
wing  box/control  surface  configurations.  A  wing 
assembly  and  a  canard  or  horizontal  tail  may  be 
attached  to  a  fuselage  (modeled  as  a  flexible  beam  or 
a  flexible  plate)  to  simulate  complete  airplane  config¬ 
urations.  The  level  of  modeling  detail  can  be 
selected  independently  for  each  section.  Therefore 
the  degree  of  detail  used  to  model  control  surfaces 
for  analysis  and  synthesis  is  not  limited,  as  is  the  case 
in  the  TSO  code. 

At  the  present  stage  of  research  structural  toplo- 
gy,  shape  and  material  properties  arc  preassigned. 
Skin  layer  fiber  orientation  arc  available  as  design 
variables,  for  skin  layer  thicknesses  (F.q.  3)  the 
coefficients  of  the  thickness  power  scries  sene  as 
design  variables.  This  guarantees  smooth  thickress 
variation  for  each  layer.  For  spar  and  rib  cap  areas 
(Fq.  4)  two  coefficients  are  used  as  design  variables 
for  each  spar  or  rib.  Concentrated  masses  at  preas¬ 
signed  locations  and  spring  constants  for  linear  and 
rotational  springs  can  also  be  treated  as  design  vari¬ 
ables. 

Aerodynamic  Design  Variables 

Wing  cross  section  or  aerodynamic  planform 
topology  arc  preassigned  here.  Performance  and 
loads  in  quasi-static  maneuvers  can  be  influenced  by 
the  jig  shape  (initial  camber)  of  the  wing  and  by- 
proper  deflection  of  leading  edge  and  trailing  edge 
control  surfaces.  The  initial  camber  of  the  wing  is 
given  by  a  scries 

a; 

(5) 

(=l 

where  the  powers  m,  and  n,  arc  identical  to  those  in 
Fq.  1  and  any  subset  of  the  coefficients  can  serve 
as  design  variables.  The  deflections  6,  of  control 
surfaces  (for  each  maneuver  point  scperately)  arc 
also  available  as  design  variables. 

Control  System  Design  Variables 

The  control  system  design  variables  of  the  sizing 
type  consist  of  the  values  of  coefficients  in  the 
numerator  and  denominator  of  control  law  transfer 


99 


factor  techniques  if  any  measured  data  are  available. 
As  to  control  system  modeling,  the  techniques  used 
here  make  it  possible  to  properly  model  real  flight 
control  or  active  flutter  suppression  ‘terns,  Fffects 
of  real  actuators  and  sensors  on  v’rol  system 
performance  arc  automatically  takci.  into  account 
through  their  given  transfer  functions. 

Desien  Variables 

Ref.  22  sets  forth  a  framework  for  multidiscipli¬ 
nary  wing  synthesis  and  describes  a  hierarchy  of 
design  variables  consisting  of  sizing  type  design  vari¬ 
ables  at  the  lowest  level,  followed  by  shape  and  then 
toplogical  type  design'  variables.  'l"his  classification 
applies  to  design  variables  spanning  struc¬ 
tures, control  and  aerodynamics.  Here,  in  addition  to 
a  balanced  approach  to  behavior  prediction  in  terms 
of  the  analysis  techniques  selected,  a  balance  is  main¬ 
tained  in  level  of  optimization  by  focusing  on  sizing 
type  design  variables  for  all  disciplines  considered; 

Structural  Design  Variables 

Figure  2  shows  an  airplane  modeled  as  an  assem¬ 
bly  of  flexible  lifting  surfaces,  Fach  lifting  surface  is 
modeled  as  an  equivalent  plate  whose  stiffness  is 
controlled  by  contributions  from  thin  cover  skins 
(fiber  composite  laminates)  and  the  internal  structure 
(spar  and  rib  caps).  Plate  sections  are  connected  to 
each  other  via  stiff  springs  (representing  hinge  stiff¬ 
ness  at  attach  points)  and  flexible  springs  (represent¬ 
ing  the  stiffness  of  actuators  and  their  backup 
structure).  F.ach  wing  section  can  be  made  of  several 
trapezoidal  parts.  Concentrated  masses  are  used  to 
model  nonstructural  items  and  balance  masses. 

The  vertical  displacement,  w,  of  each  wing 
section  is  approximated  by  a  Ritz  polynomial  series 
of  the  form 

A, 
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"here  x  and  y  are  chordwisc  and  span  wise  coordi¬ 
nates  respectively,  m,  and  n,  are  powers  reflecting  the 
«>Pv  of  polynomial  series  used.  It  can  be  a  complete 
polynomial  in  x  and  y  or  a  product  of  polynomials 
m  x  and  y  (Ref.  23).  The  depth  of  a  wing  section  is 
gi'en  by  a  polynomial 
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''here  the  //,  are  preassigned  parameters.  Thickness 
di'tnhution  of  a  typical  skin  layer  is  represented  by 
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Rib  and  spar  cap  areas  are  allowed  to  vary  linearly 
over  their  length  rj 

Mn)=  Ao  +  A,  t]  (4) 

The  present  equivalent  plate  modeling  capability25 
makes  it  possible  to  efficiently  analyze  combined 
wing  box/control  surface  configurations.  A  wing 
assembly  and  a  canard  or  horizontal  tail  may  be 
attached  to  a  fuselage  (modeled  as  a  flexible  beam  or 
a  flexible  plate)  to  simulate  complete  airplane  config¬ 
urations.  The  level  of  modeling  detail  can  be 
selected  independently  for  each  section.  Therefore 
the  degree  of  detail  used  to  model  control  surfaces 
for  analysis  and  synthesis  is  not  limited,  as  is  the  case 
in  the  TSO  code. 

At  the  present  stage  of  research  structural  toplo- 
gy,  shape  and  material  properties  arc  preassigned. 
Skin  layer  fiber  orientation  arc  available  as  design 
variables.  For  skin  layer  thicknesses  (F.q.  3)  the 
coefficients  of  the  thickness  power  series  serve  as 
design  variables.  This  guarantees  smooth  thickrcss 
variation  for  each  layer.  For  spar  and  rib  cap  areas 
(Fq.  4)  two  coefficients  are  used  as  design  variables 
for  each  spar  or  rib.  Concentrated  masses  at  preas¬ 
signed  locations  and  spring  constants  for  linear  and 
rotational  springs  can  also  be  treated  as  design  vari¬ 
ables. 


Wing  cross  section  or  aerodynamic  planform 
topology  arc  preassigned  here.  Performance  and 
loads  in  quasi-static  maneuvers  can  be  influenced  by 
the  jig  shape  (initial  camber)  of  the  wing  and  by 
proper  deflection  of  leading  edge  and  trailing  edge 
control  surfaces.  The  initial  camber  of  the  wing  is 
given  by  a  scries 

a; 

(5) 

i=l 

where  the  powers  m,  and  n,  arc  identical  to  those  in 
Fq.  1  and  any  subset  of  the  coefficients  can  serve 
as  design  variables.  The  deflections  8,  of  control 
surfaces  (for  each  maneuver  point  seperately)  arc 
also  available  as  design  variables. 

Control  System  Design  Variables 

The  control  system  design  variables  of  the  sizing 
type  consist  of  the  values  of  coefficients  in  the 
numerator  and  denominator  of  control  law  transfer 


Aerodynamic  Design  Variables 
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analyses  (e.g.  3).  Finite  dimensional  state  space 
unsteady  aerodynamic  approximations  can  be  gener¬ 
ated  for  the  full  order  aerodynamic  matrices  at  the 
start  of  optimization,  or  directly  for  modally  reduced 
aerodynamic  matrices  during  optimization  whenever 
the  modal  basis  is  changed22 . 

Gust  Response  Analysis 

'Idie  RMS  (root  mean  square)  values  of  control 
surface  deflections  and  rates  as  well  as  RMS  values 
of  selected  sensor  measurements  due  to  continuous 
atmospheric  turbulence  arc  calculated  for  different 
(light  conditions.  Both  Dryden  and  rational  approxi¬ 
mations  for  the  Von-Karman  turbulence  spectra  are 
implemented.  The  relevant  quantities  are  RMS 
values  of  control  surface  deflections  { qc },  rates  {qc) 
and  sensor  measurements  {y5£}  . 

The  state  space  equations  (P.q.  8)  are  .transformed 
into  standard  form 

s  {.m  -  ih  w,  +  ih  hu)  (9) 

Since  only  {<7,}..f{«7,}  and  {ySB)  arc  considered22  the 
output,  yt,  is  given  by 

yk={ck)T{x)  GO) 

where  {C4}  is  either  constant  or  a  function  of  control 
system  design  variables. 

The  state  covariance  matrix  is  a  solution  of  a 
Lyapunov's  matrix  equation”  in  the  form 

[A-][X]+[X)[Af  =-{D)[Qw-\{B)T  (11) 

where  [Q.]  is  the  intensity  matrix  of  the  gaussian 
white  noise.  u\  The  Ucsscnbcrg-Schur  method40  is 
used  to  solve  Gq.  1 1. 

Induced  Drag  Analysis 

Induced  drag  is  calculated  for  the  clastic  lift 
distribution  during  maneuvers.  Drag  values  assum¬ 
ing  cither  full  leading  edge  suction(  fully  attached 
flow)  or  no  leading  edge  suction  (separated  flow  at 
the"  leading  edge)  or  a  combination  of  these41  arc 
available.  The  induced  drag  D,  can  be  expresssed  in 
quadratic  form 

D,  =  ±pU2Sfi)TlADm  (12) 

where  ~pU>  is  the  dynamic  pressure  and  S  is  a 
reference  area  and  [ AD ]  is  an  aerodynamic  matrix. 


When  the  jig  shape  generalized  displacements  {<?*} 
and  control  surface  rotations  (<5)  arc  small,  then 
(wjthin  linearized  lifting  surface  theory)  this  matrix  is 
fixed  and  docs  not  depend  on  either  structural  or 
aerodynamic  design  variables42 .  It  is  generated  once 
for  the  set  of  -Rif/,  polynomials  (Hq.  1),  control 
surface  rotations  and  down  wash  due  to  pitch  rate 
degrees  of  freedom  in  symmetric  pull-up  maneuvers. 
q  is  the  following  vector 

-  K] 

{?}  =  +  <  Q  >  (13) 

lo)  (0  J 

where  <?  and  q°  arc  definded  in  Gqs.  1  and  5.  'Ihc 
vector  8  contains  all  control  surface  rotations  and  8 
is  the  pitch  rate.  Thus,  D,  depends  on  the  structural 
and  acrody  namic  design  variables  through  the  gener¬ 
alized  deflections  {<?}  in  maneuver,  the  jig  shape 
generalized  coefficients  (qn)  and  control  surface 
rotations  {8}  . 

Behavior  Functions 

The  following  behavior  functions  can  be  evalu¬ 
ated  with  the  present  capability:  clastic  displace¬ 
ments;  clastic  twist;  spar  rib  cap  stresses;  skin 
combined  stress  failure  criteria43;  natural  frequencies; 
real  and  imaginary  parts  of  acroscrvoclastic  poles; 
RMS  of  random  control  surface  deflections  and  rates 
due  to  gust;  RMS  of  sensor  measurements  in  gust; 
total  mass;  lift  and  drag  coefficients;  control  surface 
deflections  and  hinge  moments  in  maneuvers;  roll 
rate  or  load  factor  in  maneuvers. 

Control  surface  effectiveness  is  not  addressed 
directly  at  this  stage.  Instead  the  synthesis  focuses 
on  sustaining  a  desired  roll  rate  or  load  factor  while 
keeping  hinge  moments,  control  surface  deflections 
and  stresses  within  allowable  bounds. 

Acroscrvoclastic  stability  is  gauaranteed  by 
providing  enough  damping  at  each  flutter  critical 
acroscrvoclastic  pole  throughout  the  flight 
envelope44.  Handling  qualities  can  be  addressed  via 
inequality  constraints  on  the  acroscrvoclastic  pole 
locations  (e.g.  short  period  root  placement)  and  pilot 
seat  acceleration  due  to  atmospheric  turbulence45-44. 
The  control  surface  deflection  needed  for  trim  and 
overall  performance  in  a  given  maneuver  and  its 
RMS  activity  due  to  gusts  can  be  combined  in  a 
single  constraint  to  avoid  saturation4’-4*. 

Individual  behavior  functions  or  their  combina¬ 
tions  can  serve  as  objective  functions.  Possible  alter¬ 
natives  are  mass,  drag(to  be  minimized),  steady  roll 
rate  or  lift  to  drag  ratio  (to  be  maximized).  RMS  of 
aileron  rotation  or  rotation  rate  due  to  turbulence  or 
a  combination  of  any  of  these. 
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Behavior  Sensitivity  Analysis 


Test  Cases 


Implicit  differentiation  of  the  analysis  equations 
is  used  here  to  derive  analytical  expressions  for  the 
derivatives  of  all  behavior  functions  with  respect  to 
all  design  variables22,49.  This  results  in  a  computer 
code  which  is  much  larger  and  more  complicated 
than  if  finite  difference  derivatives  were  us'  l. 
However.  analytical  sensitivities  are  freefrom  numci 
ical  problems  associated  with  finite  difference  step 
si/.c  selection  and  are  attractive  in  the  context  of 
multidisciplinary  synthesis  because  of  superior 
computational  efficiency. 

Approach  to  Integrated  Optimization. 

Once  the  preassigned  parameters,  design  vari¬ 
ables,  failure  modes,  load  conditions  and  objective 
function  arc  selected,  the  integrated  optimization 
problem  can  be  cast  as  a  nonlinear  programming 
problem  of  the  form 

min,*)  F({ X}  )  (14) 


(S({.V})}  <  {0} 

{XL)<{.\-)<(XL) 

where  F  is  the  objective  function,  {/V}  is  the  vector 
of  design  variables,  {g}  is  a  vector  of  inequality 
constraints  and  (AT),  {Ax}  are  vectors  of  design  vari¬ 
able  lower  and  upper  bounds. 

The  nonlinear  programming  approach  combined 
with  approximation  concepts  (NLP  AC  approach) 
has  proven  to  be  an  effective  method  for  solving 
structural  synthesis  problems12,15  and  here  it  is 
adapted  to  the  multidisciplinary  design  optimization 
task.  In  this  method  relatively  few  detailed  analyses 
are  carried  out  during  optimization.  Each  analysis 
and  the  associated  behavior  sensitivity  analysis  serve 
as  a  basis  for  constructing  approximations  to  the 
objective  and  constraint  functions  in  terms  of  the 
design  variables.  Thus,  a  series  of  explicit  approxi¬ 
mate  optimization  problems  is  solved  converging  to 
an  optimal  design. 

The  main  advantage  of  this  approach  is  in  its 
generality.  No  apriori  assumptions  have  to  be  made 
about  the  set  of  active  constraints  at  the  optimum. 
Given  an  initial  design,  a  local  optimum  is  sought 
using  mathematical  programming  techniques.  Thus 
it  is  especially  suitable  for  multidisciplinary  optimiza¬ 
tion,  where  the  problem  is  large  and  complicated  and 
past  experience  does  not  provide  much  intuitive 
guidance.  However,  for  this  approach  to  be  practical 
it  is  crucial  to  avoid  too  many  detailed  analyses. 
Efficient  analysis  sensitivity  calculations  and  robust, 
explicit  approximations  are  essential  in  this  context. 
The  CON.M1N-'0  code  is  used  here  for  constrained 
function  minimization. 


The  accuracy ,  power  and  computational  efficien¬ 
cy  of  the  present  capability  are  discussed  in  Ref.  22 
using  a  mathematical  model  of  a  light  weight  fighter 
similar  to  the  YI  -16.  I'or  the  present  optimization 
studies  a  mathematical  model  of  a:  small  remold) 
piloted  vehicle  is  used.  Its  planform  geometry  is 
shown  in  Fig.  3.  A  6.8  aspect  ratio,  bicomex  10°o 
t/c  wing  is  actively  controlcd  by  a  small  control 
surface  located  at  about  80%  semi-span  towards  the 
tip.  The  control  surface  chord  is  20%  of  the  local 
wing  chord,  and  it  is  driven  by  an  actuator  whose 
transfer  function  is  preassigned 

fti  _  1. 7744728a:  107  (15) 

•V  {(s+  180Xr2+251r+3142)} 


The  wing  control  surfaces  are  only  used  for  active 
flutter  control.  The  elevators  are  used  for  rigid  body 
pitch  and;  roll.  The  elevator  actuator  transfer  func¬ 
tion  is 


<?<2  _  20 
<52  (r+20) 


(16) 


q„  and  5,  arc  the  actuator  actual  deflection  and  its 
command,  respectively. 


The  RPV  structure  is  modeled  as  an  assembly  of 
four  equivalent  plates.  A  flexible  wing  is  attached  to 
a  rigid  fuselage  and  rigid  control  surfaces.  The  wing 
is  divided  into  three  trapezoids.  The  two  trailing 
edge  extensions,  to  the  left  and  right  of  the  control 
surface,  arc  assumed  fixed.  The  main  wing  box  struc¬ 
ture,  extending  from  root  to  tip  spanwise  and  to 
80%  chordwisc,  is  the  structure  to  be  synthesized, 
and  alternative  designs  can  used  for  the  test  exam¬ 
ples.  The  weight  of  fuselage,  control  surfaces  and 
non-structural  wing  mass  is  308  Kgm  for  a  half 
airplane.  A  Drydcn  gust  model  with  a  scale  length 
of  518.16.mctcrs  and  a  vertical  gust  RMS  velocity  of 
1.06  m/scc  is  used. 


Following  Refs.  51,52  an  accelerometer  is  placed 
on  the  wing  strip  containing  the  control  surface.  It 
is  located  in  the  middlc(spanwisc)  and  0:65  chord 
point  of  the  strip.  Its  measurement  serves  as  an  input 
to  a  control  law  which,  in  turn,  generates  an  input 
command,  S,  to  the  actuator  of  the  wing  control 
surface. 


The  set  of  load  conditions  for  wing  stress  calcu¬ 
lations  consists  of  three  2g  symmetric  pull-ups  at  sea 
level,  10,000  feet  and  20,000  feet.  In  the  'maneuver 
load'  calculations  (Eq.  6)  the  airplane  is  trimmed 
using  the  elevator.  All  stress  constraints  reflect  a  1.5 
safety  factor.  Flutter,  gust  and  acroscrvoclastic 
stability  calculations,  though,  are  carried  out  at  sea 
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level,  Mach  0.9  for  the  cantilevered  wing  in  the 
examples  described  herein.  1'his  is  dbne  inten¬ 
tionally  in  order  to  first  examine  flutter  suppression/ 
structural  optimization  using  a  realistic  but  simple 
example,  leaving  flight  mechanics  issues  for  future 
studies. 

It  should  be  reemphasized  that  the  simplifying 
assumptions  above  are  made  only  for  illustrative 
.purposes  and  to  facilitate  physical  interpretation  of 
the  design  optimization  results.  The  present  capabil¬ 
ity  can  handle  airplane  models  where  multiple  equiv¬ 
alent  plate  elements  are  synthesized  subject  to  several 
pull-up  or  rolling  maneuvers.  Control  systems  can 
contain  many  control  elements  and  control  laws,  and 
aeroservoclastic  gust  response  analyses  can  be  carried 
out  for  symmetric'  anti-symmetric  free  free  motion 
in  several  flight  conditions.  IV sing  the -data  manager 
of  Ref.  53,  the  modeling  detail  and  model  size  are 
only  limited  by  available  computer  memory.  The 
(‘IT  time  limits  will  determine  the  number  of  load 
conditions  and  flight  conditions  for 
aeroservoclastic  gust  response  analysis.  However,  as 
shown  later,  quite  complex  problems  can  be  handled 
with  reasonable  computer  resources. 

The  control  law  used  for  this  study  is  the  Localized 
Damping  Type  Transfer  Function  (LDTIT) 
described  in  Ref.  54.  This  low  order  control  law 
provides  damping  locally"  in  the  range  of  frequen¬ 
cies  where  damping  is  needed.  Its  form  is 

'S  =  TT~T - 7  *‘'0.6Se  (17) 

(j  +  bc  s  +  cc) 


y.v.  =  "o.65f 

R  the  accelerometer  measurement  and  ac,  b(,  cc  are 
control  system  design  variables.  The  denominator 
coefficients  can  be  associated  with  equivalent  damp¬ 
ing  and  natural  frequency  a>c  of  the  control  law 

r.-=w?  (18) 

(19) 

Hius,  c,  and  b,  determine  the  center  frequency  and 
gain  peak  width  of  the  control  law  transfer  function 
while  ac  determines  the  effective  gain. 

lhc  preassigned  accelerometer  transfer  function  is 


u* - 0 o 

(s2  +  376.8.$  +  3142) 

The  I.DTTF  control  law  is  used  here  without 
compensation  for  sensor  and  actuator  transfer 
functions”. 


Results 

For  a  first  sequence  of  numerical  studies  a  wing 
tip  pod  is  added  to  the  yving  simulated  by  two  2.5 
Kgm  masses  at  the  forward  and  aft  points  of  the  tip. 
The  wing  box  construction  consists  oLall  aluminum 
cover  skins  and  there  are  no  spars  or  ribs  in  order  to 
simplify  the  model  and  introduce  as  few  structural 
design  variables  as  possible.  All  stress  earning 
capacity  is  thus  confined  to  the  skins,  which  are  held 
together  mathematically  by  the  plate  assumptions 
and  in  practice  by  an  array  of  minimum  gage 
spars;ribs  whose  stiffness  can  be  neglected.  The 
Von-Mises  yield  criteria  is  used.  ’Ihe  skin  thickness 
distribution  is  a  nine  term  polynomial  in  x  and  y, 
whose  terms  are  formed  from  of  the  polynomial 
product  (l,ar,  jr!)  (l.j;,^1)- .  There  are  thus  nine 
structural  design  variables. 

Figure  4  shows  skin  mass  convergence  histories 
for  three  synthesis  cases,  all  starting  with  a  1  mm 
uniform  skin  distribution.  In  the  first  case,  mass  is 
minimized  subject  to  stress  and  minimum  gage 
constraints  only.  Minimum  skin*  thickness  is  0.3S1 
mm  (.015  inch),  and  the  Von  Mises  equivalent  stress 
and  minimum  gage  are  constrained  at  25  points  on 
the  wing  box  (5  chordwise  x  5  spanwisc).  This  is  the 
"stress  design"  without  flutter  constraints.  In  this  case 
the  skin  mass  is  reduced  in  1 1  full 
analyses  approximate  problem  optimization  cycles 
from  4.486  to  1.743  Kgm.  The  stress  design  is  aeroe- 
lastically  unstable.  It  flutters  at  sea  level,  Mach  0.9. 
And  thus,  a  second  synthesis  is  carried  out.  The 
same  nine  structural  design  variables  are  used  and 
the  same  stress  and  minimum  gage  limitations  are 
imposed,  however  dynamic  acroelastic  stability- 
constraints  are  now  added  to  the  set  of  requirements 
that  must  be  satisfied. 

It  is  required  that  at  sea  level,  Mach  0.9  there 
should  be  at  least  4.5%  equivalent  viscous  damping 
in  the  two  lowest  frequency  poles,  and  at  least  1.5% 
damping  in  the  next  three,  (aerodynamic  poles 
which  have  very  large  damping  ratios  are  ignored 
when  the  closed  loop  system  poles  arc  ordered  by- 
damped  frequency).  In  seven  synthesis  cycles  the 
optimization  process  reduces  skin  mass  to  3.094 
Kgm.  As  expected,  a  stiffer  and  heavier  wing  is  need¬ 
ed  to  prevent  flutter.  The  constraint  which  drives  this 
design  is  a  damping  constraint  associated  with  a  flut¬ 
ter  pole  at  14.4  Hz. 

We  next  address  an  important  question:  How- 
effective  can  an  active  control  system  be  in  further 
reducing  the  minimum  structural  weight  needed  to 
satisfy  stress  and  flutter  constraints?  Three  control 
system  design  variables  are  now  added  to  the  nine 
structural  design  variables.  Wing  skin  mass  minimi¬ 
zation  is  carried  out  subject  to  stress,  minimum  gage 
and  dynamic  acroelastic  stability  constraints  starting 


from  a  I  mm  thick  uniform  skin  and  a  control  law 
of  the  form: 

b  _  2000 

ysF-  (s2  +  40.r+  10,000) 

As  shown  in  Fig.  4,  the  skinmass  for  this  case  is . 
reduced  to  1.838  Kgm.  Examination  of  the  optimiza¬ 
tion  results  reveals  that  the  driving  constraint  is  again 
the  damping  in  a  pole.  Its  frequency  is  17.03  Hz,  and 
the  final  damping  ratio  is  .067,  which  implies  that 
with  tighter  convergence  criteria  for  terminating  the 
optimization,  additional  weight  savings  could  be 
achieved  (1%  diminishing  return  on  three  consec¬ 
utive  approximate  problem  optimizations  is  used  as 
a  convergence  criterion).  In  any  event,  the  weight  in 
the  third  case  is  brought  back  almost  to  the  level  of 
the  stress  design  weight.  Convergence  is  slower.  It 
took  27  Cycles  and  27  CPU  minutes  on  the  UC1.A 
IBM  3090-4.  Also,  while  rapid  convergence  was 
achieved  for  the  cases  with  only  stnictural  design 
variables  with  move  limits  of  40%.  it  was  necessary 
to  use  10%  move  limits  when  control  system  design 
variables  were  added  in  order  to  protect  the  accuracy 
of  system  pole  approximations.  This  explains  the 
slower  convergence  of  the  third  run. 

Minimum  mass  synthesis  of  the  wing  with  struc¬ 
tural  and  control  system  design  variables,  subject  to 
stress,  gage  and  acroscrvodastic  stability  constraints 
was  tried  again.  This  time  the  initial  structural  design 
is  the  unstable  stress  design  of  the  first  case.  Initial 
values  for  control  system  design  variables  arc  now 
different.  The  optimization  is  started  with 
1 500, '(a5  +  40a  +  4000).  Figure  5  shows  the  two  skin 
mass  convergence  histories  for  the  design  starting 
from  linm  uniform  skin  and  the  one  starting  with 
the  stress  design.  Starting  with  the  1.743  Kgm  skin 
of  the  stress  design,  the  mass  minimization  progress¬ 
es  by  first  adding  mass  to  stiffen  the  wing  followed 
by  manipulation  of  the  control  system  design  vari¬ 
ables  to  stabilize  the  wing  with  the  smallest  weight 
penalty  possible.  In  fact,  the  final  skin  weighs  1.745 
Kgm,  practically  the  same  as  if  there  were  no  flutter 
constraints  at  all. 

Final  skin  thickness  distributions  for  the  stress, 
stress  +  flutter  and  stress  +  actively  controled  flutter 
designs  are  shown  in  Fig.  6. 

The  final  control  laws  for  the  two  cases  shown  in 
Fig.  5  are  151 7.7t(s2  +  5I.4/  + 15072.9)  and 
l456.7/(sJ+ 40.2r+ 15310.4)  for  the  Imm  initial 
design  and  stress  initial  design  respectively.  It  is 
interesting  to  note  that  the  numerator  terms  (effec¬ 
tive  gains)  converge  from  different  starting  points  to 
values  that  are  within  5%  from  each  other.  The 
constant  denominator  terms  in  both  final  control 
laws  are  almost  the  same  indicating  that  active 
damping  is  introduced  at  a  band  of  frequencies 


around  19.6  Hz.  This  is  intuitively  rewarding  since 
the  flutter  mechanism  seems  to  involve  a  wing  bend¬ 
ing  root  at  16-17  Hz  and  a  second  root  at  25-26  Hz. 
The  control  law  localized  damping  action  is  thus 
tuned  so  as  to  be  approximately  in  the  middle  of  this 
band.  The  width  of  the  frequency  band  is  controled 
by  the  equivalent  control  law  damping  parameter  C. 
(Fq.  19)  which  is  21%  and  16%  for  the  two  cases  in 
Fig.  5. 

The  results  described  so  far  indicate  that,  as  long 
as  controllability  and  observability  are  guaranteed,  an 
active  control  system  of  unlimited  power  can  stabi¬ 
lize  the  wing  and  avoid  essentially  all  structural 
weight  penalty  that  would  have  been  needed  to 
adiieve  this  in  a  passive  manner,  liven  when  gust 
dynamic  stresses  become  critical  in  the  stress  design  ( 
in  addition  to  the  quasi  static  stresses  included  here) 
it  is  reasonable  to  believe  that  a  powerful  control 
system  can  save  a  substantial  amount  of  structural 
weight. 

The  next  objective  is  to  study  how  a  limited 
power  control  system  effects  an  integrated  JcMgn  and 
how  structural  weight  and  control  effort  interact  in 
the  course  of  integrated  optimization.  Gust  re'ponw 
constraints  arc  now  added  to  the  previous  set  of 
constraints.  The  RMS  values  of  control  surface 
rotation  q,  and  rotation  rate  q,  serve  as  measures  of 
control  system  effort  and  limitations-'1 

At  the  minimum  mass  control  augmented  stnic- 
tural  design  of  l  ig.  5  (starting  with  the  stress  design) 
these  RMS  values  were  0.35  degrees  anil  36  6 
deg  sec.  In  two  additional  cases  these  RMS  values 
were  constrained  placing  limitations  of  varying  sever¬ 
ity  on  the  control  system.  Figure  7  depicts  three 
mass  convergence  iteration  histories  all  starling  with 
the  stress  desist  :  a)  no  bounds  on  the  RMS  </  and 
q,\  b)  RMS  ~q,  <  0.2:,  RMSq,  <  21  7  see:  c)  RMS 
q(  <  0.1s.  R.MStf,  <  10.57  sec.  Move  limits  of  10% 
were  used  for  all  three  cases.  Nine  stnictural  and  3 
control  system  design  variables  are  used  simultane¬ 
ously. 

As  Fig.  7  shows,  convergence  within  12  cycles  is 
achieved  when  the  gust  response  constraints  are 
added  to  the  stress,  gage  and  stability  constraints. 
When  control  surface  activity  is  more  restricted,  the 
final  skin  weight  is  larger.  Thus  limited  control 
system  resources  are  traded  off  against  structural 
resources  in  a  quest  for  a  balanced  multidisciplinary 
optimum  design.  This  interaction  takes  place 
dynamically  as  the  synthesis  progresses  and  is  hard 
or  impossible  to  capture  in  sequential  parametric 
studies. 

Figure  8  adds  to  our  understanding  of  interdisci¬ 
plinary  interactions  by  following  the  skin  mass  histo¬ 
ry  (normalized  with  respect  to  the  stress  design 
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from  a  I  mm  thick  uniform  skin  and  a  control  law 
of  the  form: 

8  _  2000 
ysF-  (j2  +  40.t  +  10,000) 

As  shown  in  Fig.  4,  the  skin  mass  for  this  case  is . 
reduced  to  1.838  Kgm.  Examination  of  the  optimiza¬ 
tion  results  reveals  that  the  driving  constraint  is  again 
the  damping  in  a  pole.  Its  frequency  is  17.03  I lz,  and 
the  final  damping  ratio  is  .067,  which  implies  that 
with  tighter  convergence  criteria  for  terminating  the 
optimization,  additional  weight  sayings  could  be 
achieved  (1%  diminishing  return  on  three  consec¬ 
utive  approximate  problem  optimizations  is  used  as- 
a  convergence  criterion).  In  any  event,  the  weight  in 
the  third  case  is  brought  back  almost  to  the  level  of 
the  stress  design  weight.  Convergence  is  slower.  It 
took  27  Cycles  and  27  CPC  minutes  on  the  L'CI.A 
IBM  3090-4.  Also,  while  rapid  convergence  was 
achieved  for  the  cases  with  only  structural  design 
variables  with  move  limits  of  40%,  it  was  necessary 
to  use  10%  move  limits  when  control  system  design 
variables  were  added  in  order  to  protect  the  accuracy 
of  system  pole  approximations.  This  explains  the 
slower  convergence  of  the  third  run. 

Minimum  mass  synthesis  of  the  wing  with  struc¬ 
tural  and  control  system  design  variables,  subject  to 
stress,  gage  and  aeroservoelastie  stability  constraints 
was  tried  again.  This  time  the  initial  structural  design 
is  the  unstable  stress  design  of  the  first  case.  Initial 
values  for  control  system  design  variables  arc  now 
different.  The  optimization  is  started  with 
)500/(j2+  40;  +  4000).  Figure  3  shows  the  two  skin 
mass  convergence  histories  for  the  design  starting 
from  1mm  uniform  skin  and  the  one  starting  with 
the  stress  design.  Starting  with  the  1.743  Kgm  skin 
of  the  stress  design,  the  mass  minimization  progress¬ 
es  by  first  adding  mass  to  stiffen  the  wing  followed 
by  manipulation  of  the  control  system  design  vari¬ 
ables  to  stabilize  the  wing  with  the  smallest  weight 
penalty  possible.  In  fact,  the  final  skin  weighs  1.745 
Kgm,  practically  the  same  as  if  there  were  no  flutter 
constraints  at  all. 

Final  skin  thickness  distributions  for  the  stress, 
stress  +  flutter  and  stress  +  actively  controled  flutter 
designs  are  shown  in  Fig.  6. 

The  final  control  laws  for  the  two  cases  shown  in 
Fig,  5  are  15I7.7/(;J+ 51.4; +  15072.9)  and 
I456.7;(;J  +  40.2;  +  15310.4)  for  the  1mm  initial 
design  and  stress  initial  design  respectively.  It  is 
interesting  to  note  that  the  numerator  terms  (effec¬ 
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values  that  are  within  5%  from  each  other.  The 
constant  denominator  terms  in  both  final  control 
laws  arc  almost  the  same  indicating  that  active 
damping  is  introduced  at  a  band  of  frequencies 


around  19.6  Hz.  TTiis  is  intuitively  rewarding  since 
the  flutter  mechanism  seems  to  involve  a  wing  bend- 
ing:root  at  16-17  Hz  and  a  second  root  at  25-26  1  lz. 
The  control  law  localized  damping  action  is  thus 
tuned  so  as  to  be  approximately  in  the  middle  of  this 
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(Eq.  19)  which  is  21%  and  16%  for  the  two  cases  in 
Fig-  5. 

The  results  described  so  far  indicate  that,  as  long 
as  controllability  and  obsen ability  are  guaranteed,  an 
active  control  system  of  unlimited  power  can  stabi¬ 
lize  the  wing  and  avoid  essentially  all  structural 
weight  penalty  that  would  have  been  needed  to 
at'liieve  this  in  a  passive  manner,  liven  when  gust 
dynamic  stresses  become  critical  in  the  stress  design  ( 
in  addition  to  the  quasi  static  stresses  included  here) 
it  is  reasonable  to  believe  that  a  powerful  control 
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weight. 

The  next  objective  is  to  study  how  a  limited 
power  control  system  effects  an  integrated  design  and 
how  structural  weight  and  control  effort  interact  in 
the  course  of  integrated  optimization.  Gust  response 
constraints  arc  now  added  to  the  previous  set  of 
constraints.  The  RMS  values  of  control  surface 
rotation  q,  and  rotation  rate  q,  serve  as  measures  of 
control  system  elTori  and  limitations51  ", 

At  the  minimum  mass  control  augmented  struc¬ 
tural  design  of  Fig.  5  (starting  with  the  stress  design) 
these  RMS  values  were  0.35  degrees  anil  36  6 
deg  sec.  In  two  additional  cases  these  RAIS  values 
were  constrained  placing  limitations  of  varying  sever¬ 
ity  on  the  control  system.  Figure  7  depicts  three 
mass  convergence  iteration  histories  all  starting  with 
the  stress  desian  :  a)  no  bounds  on  the  RMS  q  and 
qc\  b)  RMS  <7,  <  0.2= ,  RMSq,  <  21  7  sec:  c)  RMS 
q(<  O.F.  R.MS<7, <  lO.o'/scc.  Move  limits  of  |fl% 
were  used  for  all  three  cases.  Nine  structural  and  3 
control  system  design  variables  are  used  simultane¬ 
ously. 

As  Fig.  7  shows,  convergence  within  12  cycles  is 
acliieved  when  the  gust  response  constraints  are 
added  to  the  stress,  g3ge  and  stability  constraints. 
When  control  surface  activity  is  more  restricted,  the 
final  skin  weight  is  larger.  Thus  limited  control 
system  resources  are  traded  off  against  structural 
resources  in  a  quest  for  a  balanced  multidisciplinary 
optimum  design.  This  interaction  takes  place 
dynamically  as  the  synthesis  progresses  and  is  hard 
or  impossible  to  capture  in  sequential  parametric 
studies. 

Figure  8  adds  to  our  understanding  of  interdisci¬ 
plinary  interactions  by  following  the  skin  mass  histo¬ 
ry  (normalized  with  respect  to  the  stress  design 
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Thus  the  design  space  spanned  three  disciplines, 
namels.  structures,  control  and  aerodynamics  and 
the  blend  of  behavior  constraints  covered  strength, 
minimum  gage,  flutter,  gust  response  and  perform¬ 
ance  (induced  drag).  The  same  RPV  model  as  used 
in  Ref.  T?  served  as  a  test;  ease. 

The  design  studies  in  Refs.  17-19  contribute  to 
the  state  of  the  art  in  design  optimization  by  devel¬ 
oping  and  gaining  experience  with  synthesis  tech¬ 
niques  for  multidisciplinary  complex  problems 
involving. a  very  rich  blend  of  constraints.  They  also 
offer  better  understanding  and  a  fresh  insight  into  the 
interactions  between  the  various  disciplines  in  wing 
design  as  well  as  numerical  results  that  may  sene  as 
a  basis  for  design  tradeoffs. 

The  purpose  of  the  present  paper  is  to  add  to  the 
growing  experience  in  acroscnoclastic  wing  opti¬ 
mization  The  design  space  for  the  numerical  exam¬ 
ples  reported  here  is  made  up  of  sizing  type 
structural  and  control  system  design  variables  (Ref. 
16)  The  present  capability  can  handle  acrodsnamic 
design  variables  as  well  but  they  arc  not  used  in  the 
present  studies.  The  RRV  (Refs.  17-19)  is  used  to 
study  optimal  designs  achieved  with  different  control 
law  structures  and  the  resulting  robustness  of  these 
control  augmented  composite  wings.  The  perform¬ 
ance  of  different  complex  eigenvalue  approximations 
is  also  examined.  7'hc  new  multidisciplinary  wing 
synthesis  capability  is  then  applied  to  the  acroscrsoe- 
ias'.ic  synthesis  of  realistic  FI 6  and  X29  type 
airplanes  in  order  to  demonstrate  its  power  and 
generality  and  assess  computational  efficiency  in 
dealing  with  more  complex  aeroelastic  and  control 
system  configurations.  Simplified  handling  qualities 
requirements  are  added  to  the  set  of  constraints. 

Analytical  Modeling  Techniques 

A  unique  integration  of  analytical  modeling  tech¬ 
niques  makes  it  possible  to  bridge  the  gap  between 
those  that  arc  over  simplified  and  the  detailed  tech¬ 
niques  that  require  >oo  much  computer  time.  In  the 
structures  area,  an  equivalent  plate  analysis  is  used 
(Refs.  20,21).  It  is  integrated  with  the  PCKFM 
(Piecewise  Continuous  Kernel  Function  Method)  for 
lifting  surface  unsteady  aerodynamics  (Ref.  22).  The 
method  of  Roger  (Ref.  23)  is  used  to  generate  finite 
dimensional  state  space  approximations  for  the 
unsteady  aerodynamic  loads.  The  integrated  aeroser- 
soelastic  ssstem  is  modeled  as  a  Linear  Time  Invari¬ 
ant  (LT1)  system.  The  control  system  is  completely 
described  b>  the  location  of  sensors  and  control 
surfaces  and  by  the  transfer  functions  of  the 
sensors.control  laws  and  actuators.  Behavior  sensi¬ 
tivity  anal>sis  is  based  on  analytical  derivatives  of  all 


behavior  functions  with  respect  to  all  design  vari¬ 
ables. 


Multidisciplinary  Synthesis  Mcthodolom 
Approach:  to  Optimization 

Following  its  success  in  structural  ssnthcsis.  the 
nonlinear  programming  approach  combined  with 
approximation  concepts  (NLP, AC)  is  used  for  the 
multidisciplinary  optimization  task  (Ref:  2-4).  In  this 
method  onh  a  small  number  of  detailed  anahsc:>  are 
carried  out  during  optimization.  Each  anah  sis  serves 
as  a  basis  for  constructing  approximations  to  the 
objective  and  constraint  functions  in  terms  of  the 
design  variables.  Then,  a  scries  of  approximate  opti¬ 
mization  problems  is  solved  converging  to  the  opti¬ 
mal  design.  For  this  approach  to  be  practical  it  is 
crucial  to  avoid  too  mans  detailed  analsses  for  func¬ 
tion  evaluation  and  deris'atisc  calculations.  This 
depends  on  making  the  approximations  accurate  yet 
simple  enough  for  efficient  solution. 

Design  Variables 

Preassigned  parameters  for  the  optimization 
include  wing  planform  and  depth  distribution,  mate¬ 
rial  properties  and  structural  layout  of  the  wing  ( 
number  of  spars  and  ribs  and  their  locations) 
Control  ssstem  structure  is  also  preassigned.  Ihus 
the  number  of  sensors  and  actuators  and  their 
locations  arc  given  along  with  the  number  of  control 
lasvs  transforming  given  combinations  of  sensor 
outputs  into  control  commands.  It  is  aho  assumed 
that  the  general  form  of  the  transfer  functions  of 
sensors  and  actuators  are  given  and  cannot  be 
changed  during  optimization. 

To  take  advantage  of  multidisciplinary  inter¬ 
actions,  the  design  space  is  opened  up  to  include 
structural  design  variables,  control  ssstem  and  acro¬ 
dsnamic  design  variables  simultaneous!)  Structural 
design  variables  then  include  polynomial  coefficients 
in  the  series  describing  skin  laser  thickness  distrib¬ 
ution  over  the  wing 

/ 

i 

Additional  structural  design  variables  include  spar  rib 
cap  areas,  concentrated  masses  and  spring  constants 
( for  the  springs  representing  stiffness  of  actuator  and 
backup  structure  connecting  control  surfaces  to  the 
wing  box  or  canard  to  the  fuselage).  Control  ssstem 
design  variables  include  polynomial  coefficients  in 
the  transfer  functions  representing  control  laws. 
Aerodynamic  design  variables  include  coefficients  in 
the  polynomial  series  for  wing  initial  (jig)  shape. 
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With  the  structural,  control  and  aerodynamic 
conficuralions  preasssigned  (following  the  hierarchy 
of  design  variables  in  Ref.  16)  the  control  augmented 
structural  synthesis  problem  formulated  in  this  work 
is  a  si/ina  problem  for  the  three  disciplines.  Thus  the 
balanced  treatment  of  these  disciplines  (controls, 
aerodynamics  and  structures)  is  also  retained  in 
formulating  the  optimization  problem. 

Objective  Functions  and  Behavior  Constraints 

The  wing  can  be  synthesized  to  minimize  mass  or 
gust  response  or  maximize  performance  with 
constraints  on  stresses,  acroservoclastic  stability, 
aircraft  performance  in  terms  of  roll  rate,  drag  or 
drag  polar  specifications  and  control  system  perform¬ 
ance  in  terms  of  activity  in  gusts  and  limits  on 
control  surface  travel  and  hinge  moment.  The 
objective  function  can  be  chosen  to  be  mass,  drag, 
RMS  value  of  any  response  to  atmospheric  turbu¬ 
lence  (to  be  minimized)  or  steady  roll  rate  or  lift  to 
drag  ratio  (to  be  maximized)  or  a  combination  of 
these.  Constraints  arc  imposed  to  meet  a  combined 
stress  criterion  for  composite  skin  layers  and  a  unidi¬ 
rectional  stress  criterion  for  spar/rib  caps.  The  acro¬ 
servoclastic  system  poles  arc  forced  to  reside  in  the 
left  hand  half  of  the  complex  plane  to  guarantee 
dynamic  stability  If  not  included  as  part  of  the 
objective  function,  the  drag,  lift,  drag,  mass  or  roll 
rate  can  be  constrained  to  ensure  acceptable 
performance  Static  load  conditions  include  sets  of 
given  loads  acting  on  the  wing  or  definition  of 
airplane  maneuvers  In  the  second  ease,  static  defor¬ 
mation  and  stresses  are  calculated  to  take  trim  and 
aeroelastic  load  redistribution  into  account. 

The  nonlinear  programming  algorithm  used  for 
constrained  function  minimization  throughout  this 
work  is  the  method  of  feasible  directions  as  imple¬ 
mented  in  the  CONM1N  code  (Refs.  25  and  26). 

The  RPV  Wine 
Description 

Optimization  studies  presented  here  deal  first 
with  a  small  RPV  similar  to  the  NASA  DAST 
research  vehicle  Its  planform  geometry  is  shown  in 
Fig.  1.  A  6.8  aspect  ratio,  biconvex  10%  t  c  wing  is 
actively  controlled  by  a  small  control  surface  located 
at  about  80%  semi -span  towards  the  tip.  The 
control  surface  chord  is  20%  of  the  local  wing  chord, 
and  it  is  driven  by  an  actuator  whose  transfer  func¬ 
tion  is  preassigned  as  follows  (Ref.  27) 


ft l  _  1.7744728* IQ7 

s\  {(.$  +  180/j2+  251  j+  3142)} 

The  wing  control  surfaces  are  used  only  for  active 
flutter  control  The  elevators  arc  used  for  rigid  body 
pitch  and  roll  The  elevator  actuator  transfer  func¬ 
tion  is 

ft 2  _  20 

b  2  (s  +  20) 

qc,  and  b,  arc  the  actuator  actual  deflection  and  its 
command,  respectively. 

An  accelerometer  is  placed  on  the  wing  strip 
containing  the  control  surface.  It  is  located  in  the 
middle(spanvvisc)  and  at  the  0.65  chord  point  of  ihe 
strip.  Its  measurement,  ysl,  serves  as  an  input  to  a 
control  law  which,  in  turn,  generates  an  input 
command,  <5,,  to  the  actuator  of  the  wing  control 
surface.  The  preassigned  accelerometer  transfer  func¬ 
tion  is 

XSE  _  3142 

''0.65c  (s2+  376.SJ+  3142) 

n0Mf  is  the  actual  vertical  acceleration  at  the  meas¬ 
urement  point. 

The  RPV  structure  is  modeled  as  an  assembly  of 
four  equivalent  plates.  A  flexible  wing  is  attached  to 
a  rigid  fuselage  and  rigid  control  surfaces.  The  main 
wing  box  structure,  extending  from  root  to  tip  span- 
wise  and  to  80%  chordwisc,  is  the  structure  to  be 
synthesized.  The  weight  of  fuselage,  control  surfaces 
and  non-structural  wing  mass  is  30S  Kgm  for  a  half 
airplane.  A  Dryden  gust  model  with  a  scale  length 
of  518.16  meters  and  a  vertical  gust  RMS  velocity  of 
1.06  nvsec  is  used. 

The  set  of  three  load  conditions  for  wing  stress 
calculations  consists  of  3g  symmetne  pull-ups  at  sea 
level,  10,000  feet  and  20,000*  feet.  In  the  'maneuver 
load"  calculations  the  airplane  is  trimmed  using  the 
elevator.  All  stress  constraints  reflect  a  1.5  safety- 
factor.  Flutter,  gust  and  acroservoclastic  stability- 
calculations,  though,  are  carried  out  at  sea  level. 
Mach  0.9  for  the  cantilevered  wing.  This  is  done 
intentionally  in  order  to  first  examine  flutter 
suppression  structural  optimization  using  a  realistic 
but  simple  example  without  flight  mechanics  inter¬ 
actions. 

The  RPV  wing  box  skins  are  made  of 
glass  epoxy  laminates  (Ref.  2$).  Tiber  directions  are 
0.  90,  +  45  and  -45  degrees  relative  to  a  line  passing 
through  the  midchord  points  of  the  wing  box.  The 
skin  is  then  modeled  as  made  of  four  unidirectional 
lamina.  The  thickness  distribution  of  each  of  these 
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I.unm.i  us  dcM-ribed  by  a  nine  term  polynomial  in  x 
,md  y  a.s  in  1  lg.  I.  Thus,  there  arc  36  structural 
design  variables. 

Structural  Designs. 

Mimimum  weight  designs  with  structural  design 
variables  only  subject  to  minimum  gage  and:  stress 
constraints  ("stress  design")  or  gage, stress  and  flutter 
constraints  ("flutter  design")  are  synthesized  first 
(lig.  2).  The  stress  design  is  unstable.  The  flutter 
constraints  arc  in  the  form  of  a  2%  lower  bound  on 
viscous  damping  in  five  modes  corresponding  to  the 
lowest  frequencies.  Move  limits  of  40%  were  used; 
and  convergence  was  achieved  within  15  optimiza¬ 
tion  cycles. 

Structure  Control-Designs. 

1  he  control .sy stem  is  now  incorporated  and  wing 
mass  is  minimized  subject  to  gage, stress  and  flutter 
constraints  while  the  design  space  includes  both 
structural  and  control  system  design  variables. 

Hollowing  its  successful  application  to  the  all 
Aluminum  wing  studied  in  Ref.  17,  the  first  control 
law  used  for  this  study  is  the  Localized  Damping 
1  y pc  'transfer  Function  (LDTTF)  described  in  Ref. 
2‘i.  This  second  order  control  law  provides  damping 
locally"  in  the  range  of  frequencies  where  damping 
is  needed.  Its  form  is 

s  _  ac 
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(s  +  b(s+  cc) 

where  is  the  accelerometer  measurement  and 
a.,  b(,  ce  are  control  system  design  variables  The 
denominator  coefficients  can  be  associated  with 
equivalent  damping  (c  and  natural  frequency  o)c  of 
the  control  law 


Thus,  cc  and  bc  determine  the  center  frequency  and 
gain  peak  width  of  the  control  law  transfer  function 
while  oc  determines  the  effective  gain. 

The  LDTTT  control  law  is  used  here  without 
compensation  for  sensor  and  actuator  transfer  func¬ 
tions  (Ref.  30).  The  active  control  system  is 
assumed  to  have  no  weight  in  the  calculations 
performed  in  this  study. 

In  the  case  of  the  all  Aluminum  wing  (Ref.  17), 
when  the  second  order  control  system  was  added  to 
the  the  problem  and  design  synthesis  started  with  the 
stress  gage  constrained  (unstable)  design,  aeroservoe- 
lastic  mass  minimization  made  it  possible  to  reduce 


the  structural  mass  to  its  stress  design  value  while  the 
active  control  system  prevented  flutter.  In  the  case 
of  the  composite  wing,  however,  when  the  design 
space  was  opened  up  to  include  the  three  control 
system  design  variables  as  well  as  the  36  structural 
design  variables,  convergence  could:  not  be  achieved 
when  starting  with  the  stress  design. 

When  synthesis  starts  with  the  feasible  flutter 
design,  convergence  is  achieved  with  10%  move 
limits,  but  only  about  60%  of  the  weight  penalty 
needed  for  flutter  prevention  is  recovered  by  the 
addition  of  the  control  system  (Fig.  3). 

Limiting  the  control  system  power  by  including 
constraints  on  the  aileron*  activity  in  atmospheric 
turbulence  yields  the  same  trend  that  was  found  for 
the  Aluminum  wing:  The  more  limited  the  control 
system  is,  the  higher  the  structural  weight  penalty 
needed  for  acroscrvoclastic  stability.  But  the  fact 
that  when  the  control  system  is  unlimited  in  power, 
it  could  not  take  care  of  the  flutter  problem  (without 
structural  penalty)  was  disturbing  and  intriguing. 
This  called  for  further  study. 

Careful  examination  of  the  iteration  histories 
starting  with  the  stress  design  revealed  complex 
eigenvalue  approximations  that  were  extremely  sensi¬ 
tive  to  design  changes.  The  accuracy  of  these 
approximations  is  evaluated  when  constraint  values 
via  a  full  analysis  at  a  new  design  point  arc 
compared  to  their  values  at  the  same  point  based  on 
approximations  constructed  from  full  analysis  (and 
sensitivity  analysis)  at  the  previous  base  design  point. 
In  the  design  synthesis  cases  that  failed  to  converge 
approximation  accuracy  for  complex  eigenvalues 
could  change  dramatically  from  iteration  to  iteration 
yielding  very  good  approximations  in  some  instances 
and  substantial  errors  in  others. 

The  Rayleigh  Quotient  Approximations  (RQA). 
(Ref.  31)  improved  approximation  accuracy  but  did 
not  solve  the  convergence  problem.  It  seems  that 
for  the  given  structure  of  the  control  system  it  is 
impossible  to  regain  all  the  weight  penalty  associated 
with  flutter  stabilization.  Indeed  a  close  look  at  the 
poles  of  the  stress  design  composite  wing  reveals  two 
flutter  mechanisms  (Fig.  4).  It  appears  that  the 
second  order  control  law  cannot  stabilize  both  simul¬ 
taneously  because  of  the  narrow  range  of  frequencies 
for  which  it  is  effective.  Indeed,  the  design  does  not 
converge  since  it  fluctuates  between  these  two  insta¬ 
bilities.  When  one  is  stabilized,  the  other  may 
become  unstable. 

The  control  law  was  subsequently  changed  to  a 
first  order  low  pass  filter  of  the  form 
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This  control  law  is  expected  to  have  a  wider 
bandwidth  and  thus  be  more  effective  in  controlling 
the  two  instabilities.  Design  iteration  histories  with 
gust  response  constraints  of  varying  severity  are 
shown  in  Fig;  5.  The  first-  order  control  law  makes 
it  possible  to  regain  most  (85%)  ofthe  flutter  weight 
penalty.  The  RMS  of  control  surface  rotation  and 
rotation  rate,  however,  arc  much  higher  compared 
with  the  second  order  control  law  (Figs.  3,5)  and  this 
is  not  surprising  given  the  wider  bandwidth  of  the 
first  order  control  law. 

A  third  control  law  was  also  studied.  It  is  a 
fourth  order  law  made  up  from  the  sum  of  two 
second  order  filters.  This  control  law  is  expected  to 
be  more  effective  than  the  single  second  order  law  in 
suppressing  the  flutter  of  the  stress  design.  This  is 
because  the  form  of  this  control  law  offers  the  free¬ 
dom  to  tunc  a  second  order-filter  to  each  of  the  two 
instabilities. 

ac  dc 

d  =  — - +  __ — -  ySE 

_{s~  +  bcs+c c)  (s  +ecs  +  fc)_ 

where  ys£  is  the  accelerometer  measurement  and 
a„  be.  C',  de,  et,  f'  arc  control  s>  stem  design  variables. 
The  denominator  coefficients  can  be  associated  with 
equivalent  damping  and  natural  frequency  we  of 
the  control  law  for  each  filter 


fc=  w2c 

bc  =  2>lcculf 

ec~  %2cu)2c 

Thus,  ce,f.  be  and  et  determine  the  center  frequencies 
and  gain  peak  widths  of  the  control  law  transfer 
function  while  ae,  dc  determine  the  effective  gains. 

The  performance  of  the  three  control  laws  when 
no  gust  constraints  arc  imposed  on  the  control 
s\stcm  is  compared  in  Fig.  6.  Indeed,  the  fourth 
order  control  law  makes  it  possible  to  stabilize  the 
wing  without  any  mass  penalty  over  the  stress 
design.  Actually,  it  even  reduces  the  mass  slightly 
with  respect  to  the  stress  design,  but  this  is  associ¬ 
ated  with  the  convergence  criteria  used.  The  stress 
design  mass  could  be  slightly  reduced  by  tightening 
the  percent  change  in  objective  function  used  in  a 
diminishing  return  comcrgcncc  criterion  (  1%  was 
used  in  the  examples  shown). 

The  initial  fourth  order  law  is 

1400./ha  +  20.J+  14000.)  + 

22000,/D2  +  100.5  +  20000.) 


The  final  fourth  order  law  is 

1421.7/(r2  +  20.65+  9596.8)  + 

21224-6/(r2+  1 20.  \s+  20752:3) 

The  double  quadratic  control  law  is  thus  tuned  to 
frequencies  of  15.6  and  23.0  cps  with  equivalent 
damping  ratios  of  10.5%  and  41.7%  in  theses  two 
frequencies,  respectively. 

Freezing  the  stress  design  of  the  composite  RPV 
wing  and  looking  for  the  "best"  control  system  to 
stabilize  it  using  a  fourth  order  control  law  leads  to 
the  design  history  in  Figure  7.  Six  control  system 
design  variables  arc  now  used  in  an  effort  to  mini¬ 
mize  the  RMS  of  aileron  rotation  due  to  atmospher¬ 
ic  gusts.  10%  move  limits  arc  used  to  comcrgc  to  a 
control  law 

2289 ./(r2+  1 7.2s  +  9209.7)  + 

50900./(52  +  305.5  +  23380.) 

yielding  a  minimum  RMS  rotation  of  1,9  degrees 
and  associated;  RMS=of  rotation  rate  of  249.  deg  sec. 
When,  during  the  optimization,  stability  is  lost,  the 
-gust  response  calculation  is  b\ passed.  This  shows  up 
as  gaps  in  the  design  history  of  Figure  7  Imposing 
more  severe  gust  response  constraints  on  the  control 
system,  makes  it  necessary  to  trade  in  some  weight 
as  shown  in  Fig.  8. 

Of  course,  it  might  be  argued  that  control  system 
power  translates  in  the  end  to  added  mass,  and  that 
for  the  tradeoff  studies  to  be  more  realistic,  we  need 
to  include  the  weight  of  the  control  system  in  the 
objective  function.  The  current  capability  can  take 
this  into  account  by  linking  the  values  of  certain 
concentrated  masses  to  the  RMS  aileron  rotation 
and  rotation  rate  needed.  This  was  not  carried  out. 
however,  in  the  examples  given  here  because  of  lack 
of  appropriate  data  that  will  make  such  a  linking 
meaningful.  In  any  case,  as  with  the  Aluminum 
wing,  the  complex  tradeoff  between  structural  weight 
and  control  system  power  (or  the  resulting  weight  of 
the  control  system)  is  evident. 

Another  issue  of  extreme  importance  in  control 
system  synthesis  is  that  of  robustness  (Refs.  30.32). 
Although  robustness  is  not  included  dirccth  in  the 
set  of  behavior  functions  in  this  stud\.  it  is  interest¬ 
ing  to  examine  robustness  of  the  control  systems 
synthesized  in  light  of  the  fact  that  control  system 
synthesis  here  is  carried  out  for  a  plant  that  is  chang¬ 
ing  during  the  synthesis  process. 

All  the  examples  up  to  this  point  involve  a  single 
input  single  output  control  system.  Robustness  of 
this  control  s>stem  can  be  studied  b>  examinimg  the 
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\u|iiN  plots  the  open  loop  system.  These  are 
•.linwnrin  Figs.  9-12. 

I  nures  9-12  show  Nyquist  plots  for  control 
.iii'.'muited  wine  designs  subject  to  gage, stress  and 
Duller  constraints  for  the  quadratic,  first  order  and 
lourili  order  control  laws,  respectively.  In  studying 
i hem  it  should  be  remembered  that  they  reflect 
umirol  laws  synthesized  for  different  final  plants 
(ume  structure).  The  most  robust  is  the  second 
order  controNaw  (with  the  highest:  penalty  in  terms 
of  wing  weight).  The  first  and  fourth  order  laws 
show  poor  robustness  (small  gain  and  phase 
margins)  as  a  result  of  the  increased  bandwidth. 
Considerable  improvement  in  gain  margins  is 
achieved  by  subjecting  the  system  to  gust  response 
constraints  (fig.  12).  This  results  in  a  higher  weight 
and  corresponds  to  higher  damping  in  the  resulting 
aeroscrvoelastic  poles. 

Before  concluding  this  section,  it  is  interesting  to 
c.samine  the  effect  of  different  complex  cigcmaluc 
approximations  on  synthesis  results.  Composite  wing 
skin  mass  histories  for  a  structural  design  subject  to 
stress,  minimum  gage  and  flutter  constraints  are 
shown  in  fig.  13.  The  hybrid  approximations  lead  to 
a  slightly  heavier  design  while  RQA  and  direct 
laylor  senes  approximations  lead  to  csscntialy  the 
same  result. 

The  P16Tvpe  Airplane  Model 

The  studies  presented  thus  far  focused  on  a  very’ 
simple  airplane  configuration  (the  RPV)  and  a 
simple  single  input  single  output  control  system.  In 
ord„r  to  demonstrate  the  power  of  the  present  capa¬ 
bility  in  synthesizing  more  complex  configurations 
with  more  complex  control  systems  two  more  realis¬ 
tic  models  are  considered. 

The  first  is  similar  to  an  F 1 6  fighter  airplane  and 
is  shown  in  Fig.  14.  A  flexible  wing/flaperon  combi¬ 
nation  is  attached  to  a  rigid  fuselage,  elevator  combi¬ 
nation.  The  wing  box  skin  thickness  distribution  is 
to  be  synthesized.  The  airplane  is  initially  statically 
stable  and  weighs  4234  Kgm  per  half  airplane  (not 
including  the  skin  mass).  Minimum  gage  and  stress 
constraints  are  imposed  on  an  array  of  5  x  5  grid 
points  over  the  skin.  The  two  maneuver  conditions 
considered  for  stress  calculations  are  a  7.33g  symmet¬ 
ric  pullup  and  a  steady  160  deg, 'sec  roll,  both  at  sea 
level,  ,\1  =  0.9.  Stress  (stress  +  gage  constraints)  and 
flutter  (stress  +  gage  +  flutter)  design  histories  are 
shown  in  Fig.  15  for  an  all  Aluminum  wing  and  a 
composite  Graphite,  Epoxy  wing.  Minimum  gage  is 
0.50Smm  (0.02  inch)  for  the  Aluminum  skin  and 
0.0127mm  (  0.005  inch)  for  each  laminate  (  0,  90, 
+  45,  -45  deg.)  in  the  Gr.Ep  skin.  Aeroelastic  stabil¬ 
ity  is  examined  for  symmetric  and  anti  symmetric 


vibrations  of  the  free  free  airplane  at  sea  level. 
M  =  0.9.  A  minimum  of  5%  damping  (  £  )  is 
required  in  the  first  four  symmetric  and  three  anti¬ 
symmetric  poles  (the  first  symmetric  pole  is  a  short 
period  pole).  Damping  of  1%  is  required  for  higher 
frequency  poles. 

As  Fig.  15  shows,  convergence  of  the  design 
process  is  achieved  within  15  full 
analysis, optimization  cycles.  When  only  structural 
design  variables  were  considered,  mmc  limits  of  40% 
w’crc  used.  The  big  penalty  in  terms  of  weight  paid 
for  meeting  the  flutter  constraints  is  clearly  evident  as 
well  as  the  weight  savings  possible  with  composite 
construction. 

A  multi  input  multi  output  (MIMO)  control 
system  is  now  added  to  the  model  (Fig.  16).  It  is 
structured  after  the  actual  F16  control  system  (Ref. 
33)  with  an  addition  for  flutter  suppression.  Angle 
of  attack,  pitch  rate  and  normal;  acceleration  arc 
measured  at  the  center  fuselage  and  used  by  control 
laws  that  take  care  of  static  stability  and  handling 
qualities  of  the  airplane  in  symmetnc  motion.  Roll 
rate  measurement  at  the  fuselage  is  used  for  the  roll 
channel  to  control  rolling  performance.  Wing  tip 
accelerometers  (at  the  tip  leading  edge)  of  each  wing 
are  used  for  flutter  suppression  in  combination  with 
the  fuselage  normal  acceleration  reading.  The  sum  of 
the  tip  accelerations  is  used  for  flutter  suppression  in 
symmetric  motion.  Their  difference  is  used  for  anti 
symmetric  vibration  stabilization. 

Assuming  perfect  sensors  and  using  the  transfer 
functions  of  the  actual  F16  actuators  (Ref.  33),  six 
control  laws  are  synthesized  simultaneously  to  ensure 
5%  damping  (£  )  in  poles  associated  with  elastic 
modes  in  symmetric  and  anti  symmetric  motion.  A 
35%  minimum  damping  constraint  and  constraints 
on  short  period  frequency  that  limit  it  to  the  range 
of  0.25  -  1.00  cps  (Ref.  10)  as  well  as  a  35%  mini¬ 
mum  damping  requirement  on  the  control 
augmented  roll  pole  are  a  simple  way  to  introduce 
handling  quality  considerations.  The  interactions 
between  a  flight  control  system  and  an  active  flutter 
suppression  system  can  thus  be  taken  into  acount  in 
the  early  design  stages. 

Figure  17  shows  a  skin  mass  design  iteration 
history  for  the  control  augmented  F 16  type  airplane 
model.  A  total  of  36  design  variables  are  used  for 
the  thickness  distribution  of  a  wing  box  skin  consist¬ 
ing  of  4  composite  laminates  and  14  design  variables 
are  used  in  the  control  system.  Gage,  stress,  flutter 
and  handling  quality  constraints  are  included.  The 
design  converges  to  the  stress  design  weight  in  1 1  full 
analysis, optimization  cycles.  The  initial  and  final 
control  systems  are  shown  in  Figs.  18  and  19.  Move 
limits  of  10%  were  used. 
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Nyquist  plots  of  the  open  loop  system.  These  are 
shown  in  Figs.  9-12. 

Figures  9-12  show  Nyquist  plots  for  control 
augmented  wing  designs  subject  to  gage, stress  and 
flutter  constraints  for  the  quadratic,  first  order  and 
fourth  order  control  laws,  respectively.  In  studying 
them  it  should  be  remembered  that  they  reflect 
control  laws  synthesized  for  different  final  plants 
(wing  structure).  The  most  robust  is  the  second 
order  control  law  (with  the  highest  penalty  in  terms 
of  wing  weight).  The  first  and  fourth  Order  laws 
show  poor  robustness  (small  gain  and  phase 
margins)  as  a  result  of  the  increased  bandwidth. 
Considerable  improvement  in  gain  margins  is 
achieved  by  subjecting  the  system  to  gust  response 
constraints  (Fig.  12).  This  results  in  a  higher  weight 
and  corresponds  to  higher  damping  in  the  resulting 
acroscrvoclastic  poles. 

Before  concluding  this  section,  it  is  interesting  to 
examine  the  effect  of  different  complex  eigenvalue 
approximations  on  synthesis  results.  Composite  wing 
skm  mass  histories  for  a  structural  design  subject  to 
stress,  minimum  gage  and  flutter  constraints  are 
shown  in  Fig.  13.  The  hybrid  approximations  lead  to 
a  slightly  heavier  design  while  RQA  and  direct 
Taylor  series  approximations  lead  to  cssentialy  the 
same  result. 

The  F16  Type  Airplane  Model 

The  studies  presented  thus  far  focused  . on  a  very 
simple  airplane  configuration  (the  RPV)  and  a 
simple  single  input  single  output  control  system.  In 
order  to  demonstrate  the  power  of  the  present  capa¬ 
bility  in  synthesizing  more  complex  configurations 
with  more  complex  control  systems  two  more  realis¬ 
tic  models  are  considered. 

The  first  is  similar  to  an  F16  fighter  airplane  and 
is- shown  in  Fig.  14.  A  flexible  wing/flaperon  combi¬ 
nation  is  attached  to  a  rigid  fuselage/elcvator  combi¬ 
nation.  The  wing  box  skin  thickness  distribution  is 
to  be  synthesized.  The  airplane  is  initially  statically 
stable  and  weighs  4234  Kgm  per  half  airplane  (not 
including  the  skin  mass).  Minimum  gage  and  stress 
constraints  are  imposed  on  an  array  of  5  x  5  grid 
points  over  the  skin.  The  two  maneuver  conditions 
considered  for  stress  calculations  are  a  7.33g  symmet¬ 
ric  pullup  and  a  steady  160  deg/sec  roll,  both  at  sea 
level,  M  =  0.9.  Stress  (stress  +  gage  constraints)  and 
flutter  (stress  +  gage  +  flutter)  design  histories  are 
shown  in  Fig.  15  for  an  all  Aluminum  wing  and  a 
composite  Graphite/Epoxy  wing.  Minimum  gage  is 
0.508mm  (0.02  inch)  for  the  Aluminum  skm  and 
0.0127mm  (  0.005  inch)  for  each  laminate  (  0,  90, 
+  45,  -45  deg.)  in  the  Gr/Ep  skin.  Aeroelastic  stabil¬ 
ity  is  examined  for  symmetric  and  anti  symmetric 


vibrations  of  the  free  free  airplane  at  sea  level, 
M  =  0.9.  A  minimum  of  5%  damping  (  (  )  is 
required  in  the  first  four  symmetric  and  three  anti- 
symmetric  poles  (the  first  symmetric  pole  is  a  short 
period  pole).  Damping  of  1  %  is  required  for  higher 
frequency  poles. 

As  Fig.  15  shows,  convergence  of  the  design 
process  is  achieved  within  15  full 
analysis/optimization  cycles.  When  only  structural 
design  variables  were  considered,  move  limits  of  40% 
were  used.  The  big  penalty  in  terms  of  weight  paid 
for  meeting  the  flutter  constraints  is  clearly  evident  as 
well  as  the  weight  savings  possible  with  composite 
construction. 

A  multi  input  multi  output  (MIMO)  control 
system  is  now  added  to  the  model  (Fig.  16).  It  is 
structured  after  the  actual  F16  control  system  (Ref. 
33)  with  an  addition  for  flutter  suppression.  Angle 
of  attack,  pitch  rate  and  normal  acceleration  arc 
measured  at  the  center  fuselage  and  used  by  control 
laws  that  take  care  of  static  stability  and  handling 
qualities  of  the  airplane  in  symmetric  motion.  Roll 
rate  measurement  at  the  fuselage  is  used  for  the  roll 
channel  to  control  rolling  performance.  Wing  tip 
accelerometers  (at  the  tip  leading  edge)  of  each  wing 
are  used  for  flutter  suppression  in  combination  with 
the  fuselage  normal  acceleration  reading.  The  sum  of 
the  tip  accelerations  is  used  for  flutter  suppression  in 
symmetric  motion.  Their  difference  is  used  for  anti 
symmetric  vibration  stabilization. 

Assuming  perfect  sensors  and  using  the  transfer 
functions  of  the  actual  F16  actuators  (Ref.  33),  six 
control  laws  are  synthesized  simultaneously  to  ensure 
5%  damping  (£  )  in  poles  associated  with  elastic 
modes  in  symmetric  and  anti  symmetric  motion.  A 
35%  minimum  damping  constraint  and  constraints 
on  short  period  frequency  that  limit  it  to  the  range 
of  0.25  -  1.00  cps  (Ref.  10)  as  well  as  a  35%  mini¬ 
mum  damping  requirement  on  the  control 
augmented  roll  pole  are  a  simple  way  to  introduce 
handling  quality  considerations.  The  interactions 
between  a  flight  control  system  and  an  active  flutter 
suppression  system  can  thus  be  taken  into  acount  in 
the  early  design  stages. 

Figure  17  shows  a  skin  mass  design  iteration 
history  for  the  control  augmented  F16  type  airplane 
model.  A  total  of  36  design  variables  are  used  for 
the  thickness  distribution  of  a  wing  box  skin  consist¬ 
ing  of  4  composite  laminates  and  14  design  variables 
are  used  in  the  control  system.  Gage,  stress,  flutter 
and  handling  quality  constraints  are  included.  The 
design  converges  to  the  stress  design  weight  in  1 1  full 
analysis/optimization  cycles.  The  initial  and  final 
control  systems  are  shown  in  Figs.  18  and  19.  Move 
limits  of  10%  were  used. 
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It  is  interesting  to  notice  the  change  in  sign  in  the 
anti  symmetric  llutter  suppression  control  law,  the 
bandwidth  increase  in  the  two  flutter  suppression 
laws  and  the  practical  elimination  of  the  symmetric 
flutter  suppression  law  in  the  final  control  sy  stem 
compared  with  the  initial  design.  Locations  of  the 
complex  poles  of  the  s>mmctric  control  augmented 
FI 6  type  model  are  shown  in  Fig.  20  (actuator  poles, 
aerodynamic  poles,  integrator  poles  arc  not  shown). 
The  damping  measure  is  defined  as 


2  2 
\  ff  +  OJ 


where  a  and  u  are  the  real  and  imaginary  parts  of 
the  complex  poles.  A  positix  e  damping  measure, 
thus,  indicates  instability. 

The  synthesis  of  the  control  augmented  F 16  type 
model  wing  involved  1080  constraints  and  50  design 
variables.  It  took  18.  minutes  of  GPL'  time  on  the 
UCLA  IBM  3090  Model  600J. 

The  X29  Type  Airplane  Model 

The  X29  type  airplane  model  was  selected  for 
study  because  of  the  potential  Body  T reedom  Flutter 
(BFF)  instability  typical  of  forward  swept  wings 
(Ref.  10).  This  instability  is  a  result  of  the  inter¬ 
action  between  a  wing  bending  mode  whose  frequen¬ 
cy  drops  as  effective  stiffness  (structural  + 
aerodynamic)  is  lost  and  the  short  period  mode. 
'I  he  strong  interaction  between  rigid  body  and  elastic 
degrees  of  freedom  presents  a  serious  challenge  to  the 
control  system  designer.  A  control  system  must  now 
be  designed  that  will  ensure  proper  handling  qualities 
and  flutter  margins  simultaneously. 

The  X29  type  model  is  shown  in  Fig.  21.  A 
canard  surface  is  used  for  symmetric  trim.  The  outer 
aileron  is  used  for  roll  control.  Tip  missiles  arc  intro¬ 
duced  so  that  first  wing  bending  frequency  is  lowered 
compared  with  the  clean  wing  with  the  intention  of 
creating  a  BFF  instability  and  providing  an  interest¬ 
ing  test  case  for  the  present  studies. 

Stresses  are  calculated  in  7.33g  symmetric  pullup 
and  a  160  deg,  sec  steady  roll  at  sea  level,  M  =  0.9. 
The  skin  of  the  wing  box  is  made  of  four  layers  of 
Graphite,  Fpoxy  material  (0,  90,  +45,  -45  deg.  with 
repcct  to  a  line  connecting  the  mid  chord  point  of 
the  root  and  tip  of  the  wing).  The  same  gage  and 
stress  constraints  are  used  as  in  the  T 16  type  model 
case.  Flutter  constraints  arc  applied  to  symmetric 
and  antisymmetric  vibrations  at  sea  level,  M  =  0.9.  A 
minimum  of  3%  damping  ( £  )  is  required  in  the  first 
three  poles  associated  with  elastic  modes  and  1%  for 
higher  frequency  poles.  Six  modes  are  used  in  the 
stability  analysis. 


Figure  22  shows  skin  mass  design  histories  for 
the  stress  (stress  +  gage)  and  flutter 
(flutter  +  stress  +  v ».ge)  designs  with  structural  design 
variables  only.  Nine  term  second  order  polynomials 
arc  used  for  the  thickness  distribution  of  each  layer  ( 
as  in  the  RPV  composite  wing  case  )  for  a  total  of 
36  design  variables.  The  optimization  is  started  with 
uniform  thicknesses  and  uses  40%  move  limits. 
Convergence  is  achieved  in  10  full 
analysis/optimization  cycles. 

Figure  23  shows  a  speed  root  locus  plot  of  the 
resulting  stress  design.  As  can  be  seen,  the  stress 
design  is  unstable  and  the  instability  is  of  a  BFF 
type.  The  frequency  of  the  first  wing  symmetric 
bending  drops  and  the  short  period  root  becomes 
unstable  as  the  speed  is  increased. 

A  control  system  identical  to  the  one  used  for  the 
r  16  studies  is  now  added  and  used  as  the  starting 
point  for  the  control  system  optimization.  Skin 
mass  of  the  control  augmented  X29  type  airplane  is 
minimized  subject  to  gage,  stress,  flutter  and  simpli¬ 
fied  handling  quality  constraints  and  the  design 
history  is  shown  in  Fig.  24.  Again,  as  with  the  FI 6, 
the  control  system  takes  care  of  the  dynamic  instabil¬ 
ities  and  guarantees  proper  handling  qualities  and 
flutter  margins  while  the  structural  mass  is  reduced 
so  as  to  take  care  of  gage  and  stress  only.  Figure  25 
shows  the  final  pole  structure  of  the  control 
augmented  X29  in  symmetric  motion  (  aerodynamic 
poles,  integrator  poles  and  actuator  poles  arc  not 
shown).  1096  constraints  and  50  design  variables 
were  included  and  the  optimization  took  58  CPU 
minutes  on  the  UCLA  IBM  Model  600J. 

Examination  of  the  initial  (Fig.  18)  and  final 
(Fig.  26)  control  systems  reveals  major  changes 
introduced  during  design  optimization  including 
changes  of  sign  and  elimination  of  some  design  vari¬ 
ables.  The  power  and  applicability  of  the  present 
technology  to  complex  configurations  is  clearly 
demonstrated. 

Conclusions 


The  present  paper  and  Refs.  14-19  that  preceded 
it  clearly  show  that  using  current  supercomputers, 
ingenious  integration  of  analysis  techniques,  analytic 
sensitivities  and  approximation  concept  based  opti¬ 
mization  methodology,  the  single  level  multidiscipli¬ 
nary  synthesis  of  realistic  actively  controlled 
composite  wings  is  both  practical  and  feasible.  One 
of  the  important  lessons  of  this  research  is  that  the 
introduction  of  control  system  design  variables  in 
addition  to  structural  design  variables  presents  a 
challenge  to  current  state  of  the  art  approximation 
techniques  used  for  approximate  optimization  prob¬ 
lem  generation.  Smaller  move  limits  are  necessary 
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Abstract 

A  digital  adaptive  controller  is  applied  to  the 
active  flutter  suppression  problem  of  a  wing  under  time 
varying  flight  conditions  in  subsonic  and  transonic  flow. 
Linear  quadratic  controller  gain  at  each  lime  step  is 
obtained  using  an  iterative  Riccati  solver.  The  digital 
adaptive  optimal  controller  is  robust  with  respect  to  the 
unknown  external  loads.  Flutter  and  divergence  insta¬ 
bilities  are  simultaneously  suppressed  using  a  trailing- 
edge  control  surface  and  displacement  sensing.  A  new 
transonic  unsteady  aerodynamic  approximau'on  metho¬ 
dology  is  developed  which  enables  one  to  carry  out  the 
rapid  calculation  required  for  transonic  aeroservoelastic 
applications.  This  approximation  is  based  on  a  combi¬ 
nation  of  unsteady  subsonic  aerodynamics  combined 
with  a  transonic  correction  procedure.  Aeroservoelastic 
transient  time  response  is  obtained  using  Roger’s 
approximation,  state  transition  matrices  and  an  iterative 
time  marching  algorithm.  The  aeroservoelastic  system 
in  the  time  domain  is  modelled  using  a  deterministic 
ARMA  model  together  with  a  parameter  estimator. 
Transonic  flutter  boundaries  of  a  wing  structure  are 
computed,  in  the  time  domain,  using  an  estimated 
aeroelasiic  system  matrix  and  are  in  good  agreement 
with  experimental  data  for  the  low  transonic  Mach 
number  range. 

Nomenclature 

General  notation 

(0  Denotes  value  at  time  i. 

(»)  Denotes  steady  state  value  at  r=». 


{xh  Denotes  {*}  at  discrete  time  k. 

Scalar  Symbols 

Ai,  (»)  The  i-th  element  of  the  vector  {/!*(“)} 

Aj,(«)  The  i-th  element  of  the  vector  {A  t(~)l 

a;  AR  coefficient 

bi  MA  coefficient 

c  Chord  length 

cj  Real  part  of  the  j-th  eigenvalue 

dj  Imaginary  part  of  the  j-th  eigenvalue 

fun  Nyquist  frequency 

J  Performance  index  for  the  digital 

adaptive  optimal  controller 

2 M  Order  of  ARMA  model 

Free  stream  Mach  number 

m  Order  of  generalized  displacement 

vector  (Ti(r)J 

N,  Number  of  discrete  frequencies  at 

which  subsonic  aerodynamic  influence 
coefficient  matrices  are  computed 

NG  Number  of  aerodynamic  lag  terms 

NM  Order  of  aeroelastic  system  matrices 

obtained  from  the  Roger’s  approxima¬ 
tion 

n  Order  of  displacement  vector  {<7(r)J 
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Free  stream  pressure 

{R.) 

Generalized  transonic  aerodynamic 

<70 

Free  stream  dynamic  pressure 

force  vector  at  steady  state 

^PMl) 

{M  0} 

Generalized  transonic  aerodynamic 
force  vector  corresponding  to  the 

n 

The  i-th  element  of  the  vector  {r(i» 

unsteady  elastic  wing  motion 

r* 

Weighting  factor  for  the  control  sur¬ 
face  deflection  angle  8i 

{MO} 

Generalized  subsonic  aerodynamic 
force  vector  corresponding  to  the 

s 

Laplace  variable 

unsteady  elastic  wing  motion 

Ta 

Time  step  for  the  iterative  time- 

{ Mr)} 

Laplace  transform  of  (MO) 

marching  algorithm 

{MO} 

Generalized  transonic  aerodynamic 

T, 

Sampling  time  for  the  parameter  esti¬ 
mator 

force  vector  corresponding  to  the 
unsteady  control  surface  motion 

The  i-th  diagonal  element  of  the 
matrix  [T4(°°)] 

(MO) 

Generalized  subsonic  aerodynamic 
force  vector  corresponding  to  the 
unsteady  control  surface  motion 

>4 

Wing  response  at  sensor 

{MO} 

7 

Ratio  of  specific  heats  (=1.4) 

Laplace  transform  of  (MO) 

hm 

Control  surface  deflection  angle 

{r(OJ 

Unsteady  component  of  the  general¬ 
ized  displacement  vector 

5(*) 

Laplace  transform  of  6(0 

lXPh 

Aeroservoclastic  state  vector  associ¬ 

h 

Forgetting  factor 

ated  with  the  estimated  system  param¬ 

0; 

The  i-th  modal  damping  of  the  wing 

eters 

structure 

{*,(0) 

Aerodynamic  state  vector  correspond¬ 

«,■ 

The  i-th  aerodynamic  lag  term 

ing  to  the  elastic  wing  deformation 

0); 

The  i-lh  undamped  natural  circular  fre¬ 

{MO} 

Structural  State  vector 

quency  of  the  wing  structure 

{MO} 

Aerodynamic  state  vector  correspond¬ 

COi 

The  i-th  damped  natural  circular  fre¬ 

ing  to  the  control  surface  deflection 

quency  of  the  wing  structure 
/ 

DTO} 

Aeroservoclastic  output  vector 

c. 

The  i-lh  modal  damping  factor  of  the 

h(0) 

Generalized  displacement  vector 

Vector  Symbols 

wing  structure 

h.) 

Generalized  displacement  vector  at 
steady  state 

{0}ji 

Parameter  vector  at  lime  t=kT, 

K(01 

The  i-lh  column  vector  of  the  matrix 
K(OJ 

{0}* 

Regression  vector  at  time  i-kT, 

MiCOJ 

Generalized  transonic  aerodynamic 
influence  coefficient  vector 

Matrix  Symbols 

corresponding  to  the  unsteady  control 
surface  motion 

{AP).{BP].[CP] 

Estimated  aeroservoelasu'c  system 
matrices 

Mi(*)} 

Generalized  subsonic  aerodynamic 

lA,l.{B(y).{BIr]. 

influence  coefficient  vector 

corresponding  to  the  unsteady  control 
surface  motion 

Aerodynamic  system  matrices 

corresponding  to  the  clastic  wing 
deformation 

(Mr)} 

Laplace  transform  of  the  vector  {/i4(r)} 

{A,].{B,].IC,} 

Structural  system  matrices 

k} 

Unit  vector  whose  i-th  element  is  1 

lAiMBciMBu]. 

(G}» 

Controller  gain  vector  at  lime  i=kTt 

(Cjl.{Dci].IDIt] 

Aerodynamic  system  matrices 

{(-  )t  Estimator  gain  vector  at  time  i=kTt 

{<7(0}  Displacement  vector 

(2(0)  Transonic  aerodynamic  force  vector 

(rtf')!  Generalized  transonic  aerodynamic 

force  vector 


corresponding  to  the  control  surface 
deflection 

M,(0]  Generalized  transonic  aerodynamic 

influence  coefficient  matrix 
corresponding  to  the  unsteady  elastic 
wing  motion 
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Damping  matrix 
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Steady  state  Riccatimatrix 
Riccaii  matrix  at  time  t=kT, 

Weighting  matrix  for  the  stale  vector 
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Transonic  correction  matrix  due  to 
elastic  wing  deformation 

Transonic  correction  matrix  due  to 
control  surface  deflection 

Eigenmatrix  of  order  nxm 
Covariance  matrix  at  time  t=kT, 
Generalized  damping  matrix 
Generalized  stiffness  matrix 
State  transition  matrix  (=eIA'17*) 

State  transition  matrix  (se1*'57*) 

State  transition  matrix  (He1*4'7') 

Time  integration  of  state  transition 


matrix  [<b, 


.1 


Time  integration  of  state  transition 

T 


matrix 


ix  [«!»,]  (e^e'AlVr‘~a'do) 


Time  integration  of  state  transition 

T, 


matrix  [$>*)  0 


1.  Introduction 


2.1  Background 

The  active  flutter  suppression  problem  of 
advanced  aircraft  configurations  is  a  multidisciplinary 
design  task  which  combines  several  disciplines  such  as 
structural  dynamics,  unsteady  aerodynamics,  control, 
and  flight  mechanics.  During  the  past  decade  active 
flutter  suppression  systems  for  die  subsonic  and  super¬ 
sonic  flight  regimes  have  been  studied  extensively. 
However,  the  transonic  flight  regime  has  received  only 
very  limited  attention.  A  considerable  amount  of  effort 


has  been  spent  on  the  development  of  active  flutter 
suppression  systems  at  NASA  Langley  Research  Center 
under  the  Drones  for  Aerodynamic  and  Structural  Test¬ 
ing  (DAST)  program.1*2-3  The  majority  of  flutter 
suppression  studies  related  to  the  DAST  wing  model 
are  based  on  the  time  invariant  subsonic  aeroservoelas- 
ticity. 

Flutter  suppression  in  the  transonic  flight  regime 
using  active  controls  presents  a  challenging  problem 
due  to  the  nonlinear  nature  of  the  transonic  aeroclastic 
problem.  Simple  proportional,  integral,  and  derivative 
(PID)  type  controllers  combined  with  detailed  unsteady 
transonic  aerodynamics  based  on  computational  fluid 
dynamics  (CFD)  have  been  applied  to  the  transonic 
flutter  suppression  problem.4  In  a  recent  research 
activity,  industry,  the  Air  Force  research  laboratories, 
and  the  NASA  have  joined  efforts  on  the  Active  Flexi¬ 
ble  Wing  (AFW)  program.5  A  low-order  real  time  digi¬ 
tal  flutter  suppression  system  for  the  AFW  model  was 
actually  built  and  tested  in  subsonic  flight  conditions.6 
For  the  transonic  regime  a  simple  PID  controller  for  the 
AFW  was  studied  using  a  computational  aerodynamics 
based  simulation  7  • 8 

In  recent  years  an  extensive  body  of  research 
aimed  at  computing  the  three-dimensional  unsteady 
transonic  aerodynamic  loads  using  computational  aero¬ 
dynamics  has  been  developed.9* 10* 1 1  Excellent  com¬ 
puter  codes  based  on  the  transonic  small  disturbance 
equation  or  the  full  potential  equation  were  developed 
and  tested  with  various  wing  geometries.12* 13- 14 
Detailed  CFD  techniques  for  the  three-dimensional 
unsteady  transonic  aerodynamics  promise  accurate 
prediction  of  the  aerodynamic  forces  for  structural 
optimization  applications.  However,  these  techniques 
are  offen  difficult  to  use  in  aeroservoelastic  studies, 
because  computing  time  can  become  excessive.  Furth¬ 
ermore,  frequently  the  aerodynamic  loads  are  obtained 
in  a  form  which  is  not  convenient  for  inclusion  in  stu¬ 
dies  aimed  at  flutter  suppression. 

1.2  Motivation 

The  inherently  nonlinear  nature  of  the  transonic 
aeroelasticity  combined  with  the  high  computational 
cost  associated  with  computational  aerodynamic  codes 
has  limited  the  number  of  aeroservoelastic  studies  deal¬ 
ing  with  this  flight  regime.4*7  *8  To  reduce  computa¬ 
tional  cost  for  the  three-dimensional  unsteady  transonic 
computations,  approximation  methods  have  been 
developed  by  a  number  of  authors.15* 16* 17  Approxima¬ 
tion  techniques  developed  to  date  were  based  on  the 
frequency  domain  approach.  The  need  for  effective 
on-line  active  flutter  suppression  systems  has  led  to  the 
recognition  that  time  domain  aerodynamics  are  needed 
for  acroservoclasticity.  A  number  of  authors  have 
developed  time  domain  aeroclastic  analyses  suitable  for 
the  subsonic  or  supersonic  flow  regime.18*19 
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‘  The  majority  of  aeroservoelastic  studies  con¬ 
ducted  to  date  have  emphasized  a  linear  time  invariant 
system  approach,  because  this  approach  has  been  found 
to  be  successful  in  the  subsonic  and  supersonic  flight 
regimes.  In  the  transonic  flight  regime,  which  is 
inherently  nonlinear,  the  validity  of  this  linear  approach 
is  questionable.  Furthermore,  under  actual  flight  condi¬ 
tions  one  can  have  variations  in  flight  speed,  thermal 
loads  due  to  aerodynamic  heating,  varying  amounts  of 
fuel  in  the  wing  structure,  and  different  external  store 
configurations.  Under  such  conditions  the  time  invari¬ 
ant  assumption  may  limit  the  validity  of  the  aeroser- 
voclastic  model  used;in  the  design  of  an  active  flutter 
suppression  system.  Digital  adaptive  control  methodol¬ 
ogy  is  an  attractive  candidate  for  real  lime  active  flutter 
suppression,  because  the  control  law  parameters  can  be 
adjusted  to  follow  the  changes  in  the  aeroclastic  sys¬ 
tem. 

Discrete  time,  digital,  adaptive  control  is  a  well 
established  methodology.  The  most  attractive  feature 
of  this  approach  is  its  ability  to  handle  mildly  lime 
varying  nonlinear  system,  such  as  encountered  in  tran¬ 
sonic  aeroelasticity.  Adaptive  control  has  been  applied 
successfully  to  the  control  of  robotic  manipulators  by  a 
number  of  authors.20,21  Application  of  an  adaptive  con¬ 
troller  to  subsonic  active  flutter  suppression  has  been 
described  by  Slater  and  Livnch22  using  the  modcl- 
rcfcrcncc  adaptive  system.  The  main  drawback  of  the 
model -reference  adaptive  system  is  that  the  adaptive 
flutter  suppression  system  did  not  work  for  a  non¬ 
minimum-phase  system  which  can  occur  in  aeroclastic 
problems.22 

The  efficient  computer  modelling  and  simulation 
studies  of  a  digital  adaptive,  real  time,  active  flutter 
suppression  system  for  lime  varying  flight  conditions  in 
transonic  flow  require  the  treatment  of  three  specific 
items.  The  first  item  required  is  an  approximate 
method  for  computing  the  three-dimensional  unsteady 
transonic  aerodynamic  loads  to  reduce  the  computa¬ 
tional  cost  of  the  aeroservoelastic  studies.  Secondly, 
the  approximate  three-dimensional  unsteady  transonic 
aerodynamic  loads  should  be  calculated  in  the  time 
domain.  This  is  an  essential  requirement  for  the  real 
time  computer  simulation  of  the  adaptive  digital  flutter 
suppression  system.  The  third  ingredient  is  the  ability 
to  change  the  free  stream  Mach  number  during  the 
computer  simulation,  so  as  to  be  able  to  test  the  adap¬ 
tive  flutter  suppression  system  under  lime  varying  flight 
conditions.  This  last  ingredient  is  quite  difficult  to 
implement  because  the  unsteady  aerodynamic  load  cal¬ 
culations  are  -usually  based  on  the  assumption  of  a  fixed 
free  steam  Mach  number.  Finally,  it  should  be 
emphasized  that  a  fast,  approximate  transonic  acroscr- 
voelastic  simulation  capability  can  play  a  useful  role  in 
preliminary  design  and  structural  optimization  of 
actively  controlled  composite  wings. 


1.3  Research  Goals 

The  primary  objective  of  this  study  is  to  study 
the  application  of  a  digital  adaptive  optimal  controller 
to  transonic  active  flutter  suppression.  This  is  achieved 
by  introducing  on-line  estimation  for  the  aeroelastic 
state  space  model  combined  with  an  iterative  Riccati 
solver  for  the  control  law  synthesis  procedure.  On-line 
estimation  of  the  aeroservoelastic  state  space  model  is 
based  on  a  deterministic  Auto-Regressive  Moving 
Average  (ARMA)  model  and  an  on-line  parameter  esti¬ 
mation  technique.  Wing  response  information  at  a 
selected  sensor  location  is  used  for  on-line  parameter 
estimation. 

The  major  advantage  of  this  technique  is  that  the 
same  flutter  suppression  system  can  be  used  for  the 
subsonic,  transonic,  or  supersonic  flight  regimes,  since 
the  adaptive  control  technique  is  based  only  on  the 
time  history  of  wing  motion.  The  on-line  estimation  of 
the  aeroelastic  state  space  model  yields  additional 
benefits.  Until  now  the  least  squares  curve  fitting 
method23  and  the  moving-block  analysis24  have  been 
used  to  identify  aeroelastic  system  damping  and  fre¬ 
quencies  in  time  domain.25,26,27  However,  because  of 
excellent  convergence  of  the  parameter  estimator  used 
in  the  present  research,  the  time  domain  flutter  boun¬ 
dary  identification  procedure  can  be  carried  out  with  a 
relatively  small  number  of  aeroelastic  transient  time 
history  calculations.  Simple  eigenanalysis  can  be 
applied  to  the  estimated  state  space  equations  to  obtain 
the  aeroclastic  system  damping  and  frequencies. 

The  second  objective  of  this  study  is  to  develop  a 
simple  methodology  for  approximating  the  three- 
dimensional  unsteady  transonic  aerodynamic  loads  in 
the  time  domain  for  the  aeroservoelastic  applications. 
The  approximation  procedure  for  obtaining  these 
unsteady  transonic  loads  is  based  on  using  relevant 
information  from  lime  domain  unsteady  subsonic  aero¬ 
dynamics.  Methods  developed  for  subsonic  aero¬ 
dynamic  load  calculations  such  as  frequency  domain 
lifting  surface  theories25,29  and  the  finite  state  approxi¬ 
mations30,31,32  arc  utilized  in  this  approximation 
method.  The  transonic  approximation  method 
developed  here  requires  only  the  computation  of  a  tran¬ 
sonic  correction  matrix  which  is  used  in  conjunction 
with  a  lime  domain  subsonic  code. 

The  well  known  transonic  bucket  effect  can  be 
captured.  Thus,  transonic  aeroelastic  compuiauons  can 
be  carried  out  with  low  computational  cost  while 
retaining  the  important  physical  transonic  acroclasuc 
characteristics.  Divergence  analysis  can  also  be  earned 
out  with  considerable  case.  Another  important  feature 
of  the  approximate  transonic  aerodynamic  load 
representation  used  in  this  study  is  the  ability  to  vary 
the  free  stream  Mach  number  during  the  simulation. 
This  enables  one  to  test  the  adaptive  optimal  controller 
under  time  varying  flight  conditions. 
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2.  Mathematical  Mode! 
for  the  Aeroelastic  System  in  Transonic  Flow 

2.1  Structural  Modelling 

The  small  perturbation  equations  of  motion  for  a 
wing  structure  have  the  following  form 

IMH*<0)  +  IC]{<?(0)  +  IKltefO)  =  {<2(0}  -  (1) 

Mass  and  stiffness  matrices.  [MJ  and  (K),  in  Eq;(l)  are 
obtained  using  the  Lifting  Surface  Structural  Analysis 
(LISSA)  computer  code.33  This  code  is  based  on  wing 
equivalent  plate  analysis.34  In  this  wing  equivalent 
plate  analysis  wing  deflections  and  mode  shapes  are 
accurately  predicted  using  a  relatively  small  number  of 
polynomial  terms.  It  was  reported34  that  wing 
equivalent  plate  analysis  is  30  and  60  limes  faster  than 
the  finite  element  method  in  static  and  free  vibration 
analysis,  respectively.  The  LISSA  code  is  particulary 
efficient  for  the  multidscipiinary  optimization  of  com¬ 
posite  wings35-36  with  a  rich  variety  of  practical  con¬ 
straints. 

Solving  the  free  vibrau'on  problem  represented  by 
Eq.(l)  using  an  orthonormal  coordinate  transformation 
(<?(0)  =  (L:){q(0)  together  with  an  assumed  medal 
damping  yields: 

WO)  +  2(C<o]|Ti(0)  +  [coz]{T](r)}  =  {*(/)}  (2) 

The  orthonormalized  equations  of  motion,  Eq.(2),  can 
be  rearranged  to  obtain  the  structural  state-differential 
equation  and  structural  measurement  equation.  The 
structural  state-differential  equation  in  the  first  order 
state  variable  form  is  given  by: 

IA\(O)=iA,HX,<r)}+[B,H/?0))  (3) 

and  the  structural  measurement  equau'on  can  be  written 
as: 

0,(0}  =  [CJUl(r)j  (4) 

2.2  Time  Domain  Unsteady  Transonic  Aerodynamics 

Since  the  steady  state  transonic  aerodynamics  of 
the  wing  is  essentially  nonlinear,  the  airfoil  shape  of 
the  wing  determines  both  the  position  and  strength  of 
the  steady  state  shock.  The  perturbed  unsteady  aero¬ 
dynamics,  due  to  wing  mouon,  is  strongly  affected  by 
the  steady  state  shock  waves.  To  obtain  a  good 
approximation  to  the  unsteady  aerodynamic  loads,  the 
correct  computation  of  the  steady  state  pressure  on  the 
wing  surface  is  a  basic  requirement.37  In  some  cases 
steady  state  pressure  distribution  from  the  wind  tunnel 
test  or  computations  based  on  CFD  codes  were  used  to 
obtain  the  nonlinear  steady  state  aerodynamics  for  the 
approximate  methods.16-17,37  The  second  assumption 
frequently  used  in  approximate  transonic  aerodynamic 
load  calculations  is  that  the  perturbed  unsteady  wing 
motion  around  the  steady  state  position  can  be  treated 


as  linear  when  sufficiently  small  wing  motions  occur.38 

In  this  study,  we  adopt  the  assumption  that  for  a 
small  dynamic  motion  ofnhe  wing,  the  unsteady  tran¬ 
sonic  aeroelastic  model  is  a  dynamically  linear 
model.38  Using  this  assumption  for  the  perturbed  wing 
mouon  the  following  vector  equation  can  be  written  for 
the  generalized  transonic  aerodynamic  force  vector 
{*«>}. 

{/?(/))  =  {/?,)  +  {/?«(/)}  +  («t(/)}  (5) 

The  vector  {/?,}  in  Eq.(5)  represents  the  generalized 
steady  state  transonic  aerodynamic  force  vector. 
{!?,(/))  and  {fl{(r)I  are  generalized  unsteady  transonic 
aerodynamic  force  vectors  due  to  elasuc  deformations 
of  the-wing  around  the  steady  st2te  position  and  control 
surface  motion,  respectively. 

The  approximate  method  for  calculating  the 
three-dimensional  unsteady  transonic  aerodynamic  loads 
employed  in  this  research  consists  of  three  steps. 

1)  Steady  suite  transonic  aerodynamic  influence 
coefficient  matrices  are  computed  using  detailed 
three-dimensional  full  potential  transonic  CFD 
code. 

2)  Transonic  correction  matrices  are  obtained  from 
the  steady  state  subsonic  and  transonic  aero¬ 
dynamic  influence  coefficient  matrices. 

3)  Three-dimensional  unsteady  transonic  aerodynamic 
loads  are  obtained  from  the  three-dimensional 
unsteady  subsonic  aerodynamic  loads  corrected  by 
the  using  transonic  correction  matrices. 

The  first  step  is  the  computation  of  the  steady 
state  transonic  aerodynamic  force  vector  {/?,}.  The 
vector  [R,  1  is  a  generalized  nonlinear  steady  state  tran¬ 
sonic  aerodynamic  force  vector.  The  full  potential 
code  developed  by  Shankar  et  al.13  is  used  to  obtain 
{/?,}  and  {tj,  J.  Here,  the  generalized  displacement  vec¬ 
tor  fa, }  of  order  m  corresponds  to  the  steady  st3te 
position  of  the  wing  structure. 

The  unsteady  components  {/?,(r)J  and  {/?5(r )}  in 
Eq.(5)  are  computed  as  follows.  A  perturbed  general¬ 
ized  displacement  vector  for  the  dynamic  mouon  of  the 
wing  can  be  defined  by 

{r(r)}H{il(0}~{Tl,}  - 

Since  (r(i) }  and  the  control  surface  deflection  6 (i)  arc 
small,  it  is  assumed  that  vectors  {/?,(/)}  and  {£*(*))  in 
Eq.(5)  can  be  calculated  using  the  transonic  aero¬ 
dynamic  influence  coefficient  matrices  given  by 

{fl,<0Js(,U0J{r(0)  (6) 

and 

{*i(0}  s  lA4(r))5(r)  (7) 

{A,(OJ  and  {Ai(r)}  arc  assumed  constant  3nd  arc  calcu¬ 
lated  by  perturbing  the  generalized  displacements  about 
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die  steady  state  reference  position  of  the  wing  and  cal¬ 
culating  the  resulting  perturbations  in  transonic  aero¬ 
dynamic  forces.  Thus, 

(8) 

an 

and 

=  (9) 

and  since  there  is  no  dependency  on  time  we  denote 
them  ( A,(~))  and  {/U(~)}.  Now,  corresponding 
influence  coefficient  matrices  due  to  small  steady  dis¬ 
placement  and  control  surface  perturbauons  in  subsonic 
flow  are  calculated  by  the  Piecewise  Continuous  Kernel 
Funcu'on  Method29  (PCKFM)  and  denoted  {A<  («)]  and 
MiWl.  A  correction  matrix  is  now  constructed  to 
transform  the  subsonic  steady  generalized  load  matrix 
to  its  transonic  small  perturbation  equivalent. 

IA,(~))  =  [T,(«)]&(~)]  (10) 


leading  to: 

lTt(oo)l  = 


Similarly,  the  i-th  element  of  the  other  transonic  correc¬ 
tion  matrix  (Tt(oo}J,  which  is  a  diagonal  matrix,  is 
assumed  to  be 


This  correction  is  applied  at  each  time  r  to  unsteady 
aerodynamic  forces  calculated  by  the  PCKFM. 

KtO]  =  (t,(~)]&<0) 

(A«(01  ~fT«(~)ltf«fr)) 


leading  ( Eqs.(6)  and  (7)  )  to 
(R,(0)  =  (T,(oo)){£f<r)} 


The  second  step  in  obtaining  the  approximate  aero¬ 
dynamic  loads  used  in  this  study  resembles  the 
quasistcady  correction  method  suggested  by  Zwaan.16 

The  third  step  required  for  the  computation  of  the 
approximate  aerodynamic  loads  consists  of  the  compu¬ 
tations  of  the  generalized  subsonic  aerodynamic  force 
vectors  {/?,(»)}  and  {^(0)  in  the  time  domain.  The 
generalized  subsonic  aerodynamic  force  vectors  (#,(x)} 
and  (£tfx)}  in  the  time  domain  are  obtained  from  the 
subsonic  aerody’namic  si2ic  space  descriptions.  Using 
the  Roger’s  approximation  to  the  Laplace  transformed 
subsonic  generalized  aerodynamic  influence  coefficient 
matrices  leads  to 

+  .  (12) 

i.l 


Since  for  subsonic  unsteady  flow 

(R.(r)l  =  (A.(r)l{r(*)}  .  (13; 

wc  substitute  Eq.(13)  into  (12)  to  get 

sc 

{ft(0)=CDo-M^C»)WOIr]*(f(r)l+BCil-4r(rtj)l.(14j 

/■I  S+\lf 

Aerodynamic  states  are  defined  as 
{*„(*>)  =  {F(r)l-4- 

slXAs))  =  -n.dHX.fr))  +  ^[I]{r(f)}  1=1.2.  •  -  •  fiG 
Thus: 

r{X,(r)}  =  [A,]{*,(r)}  +  [B„Jr{FCr))  .  (15) 

where  {£(r)}T  =  [{£,( s))TlZ&))T  •  ■  •  &**(*) )7  J- 

Transforming  Eqs.(14)  and  (15)  back  to  the  time 
domain  leads  to 

{&(r)}  =  {C,]{X,(x)}  +  lD<vl{r(r)J  +  IDSrJfF</)J  .  (16) 

where  (C,]  =  [[C.HC,]  -  -  -  [C*c]). 
and  to 

{X(r»=(A,JPf,(r)}-i-{B:,J(r(0)  .  (17) 

Similarly,  for  the  control  surface  subsonic  unsteady- 
aerodynamic  forces: 

{*t(0}  =  ICdPfrfr}}  +  (D«}6(r)  +  {DI(J6{r)  (IS; 

{*i(0)  =  (AiH^(0)  +  {Bit)5(r)  (19; 

To  complete  this  section  a  brief  discussion  of  the 
limitations  of  the  present  approximations  to  the 
unsteady  transonic  aerodynamic  loads  is  presented. 
Three  possible  cases  can  be  encountered  during  the 
simulations.  The  first  case  corresponds  to  a  subsonic 
free  stream  Mach  number,  which  is  sufficiently  low  so 
that  the  local  shock  has  not  developed  yet  over  the 
wing  surface.  For  this  case  the  subsonic  assumption 
for  representing  the  unsteady  aerodynamic  effects  is 
appropriate.  The  second  flight  condition  corresponds  to 
the  development  of  a  local  shock  over  the  wing  sur¬ 
face,  however  its  posiu'on  is  still  far  upstream  of  the 
trailing  edge.  In  this  case  the  wing  trailing  edge  is  still 
in  subsonic  flow.  However,  due  to  the  local  supersonic 
subregion  on  the  wing,  the  unsteady  subsonic  assump¬ 
tion  is  violated  in  this,  fairly  small  region.  Tne  third 
condition  corresponds  to  a  strong  shock  which  is  in  the 
vicinity  of  the  trailing  edge.  Under  such  conditions 
wing  motions  will  produce  chordwise  motion  and  oscil¬ 
lations  of  the  local  shock  wave.  When  the  local  shock 
wave  moves  aft  of  the  trailing  edge,  the  wing  trailing 
edge  will  be  in  the  local  supersonic  flow-.  Obviously 
this  condition  will  represent  the  most  severe  violation 
of  the  subsonic  assumption,  used  to  generate  the 
unsteady  transonic  aerodynamic  loads,  employed  in  this 
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study. 

Anoti-;r  limitation  of  the  present  approximate 
method  is  due  to  the  quasisteady  correction  for  tran¬ 
sonic  effects.  Due  to  this  limitation  the  approximate 
method  gives  more  accurate  results  in  the  low  reduced 
frequency  range. 

2.3  Iterative  Time  Marching  Structure/Aero¬ 
dynamic  Integration 

Structural  and  aerodynamic  state-differential 
equations  are  integrated  simultaneously  using  an  itera¬ 
tive  time-marching  algorithm  and  state  transition 
matrices.  It  is  assumed  that  vectors  {/?(/)},  {r (/)},  and 
(r(O)  are  linear  in  time  interval  kTait<(.k+\)Tt.  In  this 
time  interval,  these  vectors  are  assumed  to  be  equal  to 
average  of  their  value  over  the  interval.  Where,  T. 
represents  the  time  step  of  the  iterative  time-marching 
algorithm.  Thus,  time  integration  of  state-differential 
equations  in  the  first  order  state  variable  form  yields 
the  following  state-difference  equations.25 

{x, }*„  =  i*,](x, }*  ■»  [e,][B,J  (*h  \{R]k "  .  (20) 

Similarly  the  aerodynamic  state-difference  equations 
can  be  written  as: 

(*,)»..  =  [*,)(X,h  +  [6,][B  (21) 

=  [*•){*»)*  +  [et](Bit)'5’  +2bl"  (22) 

Structural  and  aerodynamic  measurement  equations. at  a 
discrete  time  A+l  lend  themselves  to  simple  representa¬ 
tions.  The  structural  measurement  equation  given  in 
Eq.(4)  can  be  written  in  the  discrete  time  as 

(n*-i  =  (C, ]{*,}*-,  (23) 

and  aerodynamic  measurement  equations  become 

(#.)*.!  =  [C,  ){*,)*„  +  lDo,]{rh„  +  [D,r]{r  )**,  (24) 

and 

=  [CJdf,)*.,  +  {DoiJS*.,  +  {Dis)6*tl  (25) 

Jn  Eq.(20)  the  structural  state  vector  (X,)J+1  is  a  func- 
tioii  of  the  structural  state  vector  ( X ,  )4  and  the  general¬ 
ized  transonic  aerodynamic  force  vectors  { R  )*  and 
{>?  )i.|.  However,  the  generalized  transonic  aero¬ 
dynamic  force  vector  {/?}»,!  is  not  available  at  the 
beginning  of  time  step  *+l.  At  time  Jfc+l,  the  general¬ 
ized  subsonic  aerodynamic  force  vectors  can  be  written 
as 

(*.  h.i  2{Rt  )4  -  {/?.  )*.,  (26) 

{Hih.i  =  2(*,h  -  {/?»}*-.  (27) 

Thus,  Eqs.(20)  through  (25)  can  be  solved  iteratively  at 
lime  step  *+l,  and  the  iterative  time-marching  algo¬ 


rithm  can  be  summarized  as  follows: 

1)  Predict  and  using  Eqs.(26)  and 

(27). 

2)  Calculate  the  generalized  transonic  aerodynamic 
force  vector  at  time  k+ 1. 

{*')».  i  =  {/?,)  +  IT*  («)){/?,  }4.,  +  D*(~)](A»)i., 

3)  The  structural  state  vector  and  aeroservoe- 

lastic  output  vector  can  be  obtained  from 
Eqs.(20)  and  (23),  re  :c.,uvely. 

4)  Determine  the  control  surface;  deflection  angle  6,., 
and  construct  generalized  displacement  vectors 
{/•}*„,  and  {r)4.|.  The  control  surface  deflection 
angle  8^)  can  be  calculated  using  informauon 
associated  with  the  controller,  which  will  be 
described  in  the  next  section. 

5)  Calculate  aerodynamic  state  vectors  [X,  )4.,  and 

{*,)*.,  from  Eqs.(21)  and  (22),  respectively,  and 
update  generalized  subsonic  aerodynamic  force 
vectors  {£,  )»„i  and  at  time  Jfc+1  usir.g 

Eqs.(24)  and  (25).  If  the  iteration  number  is  larger 
than  the  assigned  value,  go  to  step  1.  Otherwise, 
return  to  step  2. 

Ii  should  be  noted  that  one  iterauon  in  the  this  iterative 
scheme  corresponds  to  the  time-marching  algorithm 
developed  by  Edwards  et  al.25  Two  consecutive  itera¬ 
tions  in  this  iterative  algorithm  correspond  to  the 
predictor-corrector  scheme  used  by  Robinson  et  al.39 
The  main  drawback  of  the  above  iterative  scheme  is 
that  the  time  consuming  control  law  design  procedure 
is  inside  the  iterative  scheme.  Therefore  for  better 
computational  efficiency  control  surface  deflection 
angle  84.j  is  not  updated  during  these  iterative  compu¬ 
tations. 

3.  Adaptive  Controller 

3.1  Input-Output  and  State-Space  Descriptions  of  the 
Aeroservoelastic  System 

In  this  study  single-input  single-output  deter¬ 

ministic  ARMA  model  with  2M  Auto-Regressive  (AR) 
and  2 M  Moving  Average  (MA)  coefficients40  is  used  to 
describe  the  input-output  relation  for  the  aeroservoelas¬ 
tic  system. 

2 u  w 

y *  +  'jEfityt-i  =  £8/ 5*.;  , 

i-i  /• i 

In  vector  form, 

y*  =  (0hT(dh  .  (28) 

where 

(0}*"=  |“0|  -Ol  -02 M  b\  i>2  &2M  j 

(Olfs  Jy*-1  yt-  2  ■  •  •  y4-2W  Si-;  S4_2  ■  •  ■  54-2JV  J* 
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In  Eq.(28)  yk  and  8*  are  the  wing  response  and  control 
surface  deflection  angle  at  disccte  time  k,  respectively. 
The  input-output  description  given  by  Eq.(28)  is 
equivalent  to  a  state-space  description  which  can  be 
written  as 


{XPh.i  =  (Ap]{Xp)*  +  {Bp}84 
y*  =  (Cp){XP}*  . 


(29) 

(30) 


where 


I  A,t 


“t3  j  1  0‘ 

■  0 

=  -a 2  0  1  • 

■  0 

(BP)=- 

b2 

-fly/.,  0  0 

-a-iM  0  0 

■  1 

•  0 

A2A/-1 

but 

(Vi  iXk±  m»-i{6h(<!>)nvh-> 

*4  “  tf(*4  +  ^hTlVh.,{6h  )  ' 

where  a  covariance  matrix  [V]4  is  defined  as 
[VJ*  =  [£/ J* [£>  )*[£/]*".  Detailed  descriptions  of  the  cal¬ 
culations  of  a  diagonal  matrix  [£>)t  and  an  upper  tri¬ 
angular  matrix  (i/]i.  whose  diagonal  elements  equal  to 
1,  are  presented  in  reference  47.  The  BUD  algorithm 
is  a  modified  version  of  the  recursive  least  squares 
algorithm  which  ensures  the  positive  definiteness  of  the 
covariance  matrix -in  the  recursive  estimation  procedure. 

The  main  advantage  of  the  parameter  estimator  is 
that  the  additional  state  estimator  is  not  required  to 
estimate  system  states.  Once  the  AR  and  MA 
coefficients  in  matrices  [Ap]  and  (Bp)  are  estimated, 
then  the  state  vector  [Xp)k  can  be  calculated  using 
Eq.(32). 


and. 

[Cp]=  [l  0  0  o]  . 

The  state  vector  (Xp  )*  is  defined 


(XPh 


A 

><*) 


with 


3.3  On-Line  Control  Law  Design 

The- unknown  aeroservoelastic  system  parameters 
are  estimated  using  an  on-line  recursive  estimation 
method.  The  next  step  is  an  on-line  control  law  design 
procedure  using  the  estimated  aeroservoelastic  system 
parameters.  The  estimated  aeroservoelastic  system 
parameters  are  treated  as  if  they  are  true;  i.e.  the  uncer¬ 
tainties  of  the  estimation  are  not  considered. 

The  adaptive  optimal  control  law  is  designed  to 
minimize  approximately  the  following  linear  quadratic 
performance  index  J. 


>'i  =  -  oiyi-i  +  £]8jt-i  +  A j(A— 1) 

A i(A )  =  -  ayyk. 1  +  b2 8*-i  +  A2(A— 1) 

(32) 

h-iM-i{k)  -  ~  tfyr-iyr-i  +  hyr- 164-1  +  hw.x(Jc~\) 

=  -  CyiJU-l  +  biM§k-\  ■ 

The  state  space  description  in  Eqs.(29)  and  (30)  is  in 
observer  form  and  is  completely  observable.41  Thus, 
each  cement  of  the  state  vector  given  in  Eq.(31)  can 
be  observed  completely. 


J  =  £({Xp)*r[Ql(XP}*  +r.84J)  . 

4*0 

where  [Q]  is  a  positive  semidefinite  symmetric  matrix 
and  rw  is  a  positive  constant.  The  following  control 
law  would  be  optimal  for  7. 20 

64  s-lGJflX,)*  .  (34) 

where  a  controller  gain  vector  {G}4  is  based  on 
(Gir 

rw  +  (Bp  )r[P)i  (Bp } 


3.2  On-Line  Parameter  Estimation 

Identification  of  structural  dynamic  properties 
using  parameter  estimation  techniques  has  been  studied 
by  many  authors.20, 21  • 42,43  Parameter  estimation  tech¬ 
niques  have  been  also  applied  to  the  stability  and  con¬ 
trol  studies  of  flight  vehicles.44,45,46  Application  of 
parameter  estimation  techniques  to  aeroservoelastic  sys¬ 
tems  is  introduced  in  this  study.  The  AR  and  MA 
coefficients  in  matrices  (Ap)  and  (Bp)  are  estimated 
using  the  Bierman’s  U-D  (BUD)  algorithm47  given  as 

(9)4  =  (6)4.,  +  {L)t(yk  -{0)/.,(6}») 


(Mi 


(V)i-i  (6)4 

>.4  +  (0)i7(V]1.1(0)4 


(33) 


and  the  Riccati  matrix  [PJ*  is  obtained  using 


(P]»=(Apf([P)4 


(P)t.,(Bp)(Bp}r(P)|., 

r^+{I3p}7'[P)i_J{Bp} 


)|AP]+(Q) 


(35) 


The  Riccati  matrix  (PJ4  will  converge  to  a  constant 
nonnegative  symmetric  matric  (P)  ,  which  satisfies  the 
discrete-time  algebraic  Riccati  equation 


(P)  =  (Ap)r((P)- 


[P)(Bp)(BP)r(P] 
r~  +  (Bp  )7(P)(Bp ) 


)(AP]  +  (Q) 


(36) 


,  if  the  aeroservoelastic  system  is  stabilizabie.  At  u'me 
k  the  Riccati  matrix  [PJ*  in  Eq.(35)  can  be  iterated  as 
follows; 


1)  Calculate  [P]4  based  on  Eq.(35). 

2)  If  the  Riccati  matrix  [P]*  is  converged,  go  to  the 
next  time  step.  Otherwise,  assign  [P]4.t  =  (PJ* 
and  return  to  step  1). 

Only  one  iteration  is  recommended  at  each  discrete 
time  step  to  save  the  computation  time  for  the  control 
law  design  procedure.40  The  adaptive  optimal  control 
law  given  in  Eq.(34)  is  thus  obtained  using  estimated 
aeroservoelastic  system  matrices  [Ap  ]  and  {Bp  J  together 
with  one  iteration  of  Eq.(35). 

3.4  Special  Features  of  the  Computer  Simulation 

Covariance  Matrix  Resetting  Technique. 
When  diagonal  elements,  which  corresponds  to  AR 
coefficients,  of  the  covariance  matrix  [V]t  in  the  param¬ 
eter  estimator  given  in  Eq.(33)  converge  to  very  small 
values  while  other  elements,  which  corresponds  to  MA 
coefficients,  remain  large,  then  additional  recursive 
computations  can  not  improve  the  convergence  of  AR 
and  MA  coefficients.  A  covariance  matrix  resetting 
technique  is  essential  in  order  to  increase  the  conver¬ 
gence  rate  and  robustness  of  the  parameter  estimator.41 
In  such  a  technique  the  covariance  matrix  is  periodi¬ 
cally  changed  to  a  diagonal  matrix  a[I].  In  this  study, 
the  average  of  the  maximum  and  minimum  diagonal 
elements  of  the  covariance  matrix  is  used 

(DlyalU 

({/]/  =  m  . 

where 

_  +  MiN(D;;(j-\)) 

a  2 

The  subscript  j  and  D„(J~\)  represent  the  resetting  time 
and  the  i-th  diagonal  element  of  the  matrix  [£>],-i, 
respectively. 

Order  of  the  ARMA  Model.  It  is  possible  to 
estimate  the  order  of  the  ARMA  model  from  input  and 
output  data  if  a  lattice  filter  is  used  as  a  parameter  esti- 
mator.21-43-48  The  BUD  algorithm  for  parameter  esti¬ 
mation  in  this  study  does  not  have  the  ability  to  esti¬ 
mate  the  order  of  the  ARMA  model.  Thus,  it  is 
assumed  that  the  order  of  the  ARMA  model,  2 M,  is 
known.  The  order  of  aeroelastic  system  matrices 
obtained  from  the  Roger’s  approximation  can  be  calcu¬ 
lated  from  the  following  equation.49 

NM  =  7m  +  mxNG  . 

Here,  m  and  NG  represent  the  order  of  the  generalized 
displacement  vector  and  the  number  of  aerodynamic 
lag  terms  in  the  Roger’s  approximation,  respectively. 
Theoretically,  the  order  of  a  deterministic  ARMA 
model,  2M ,  should  be  equal  to  NM.  However,  the  full 
state- feedback  control  law  and  the  small  sampling  time 


T,  for  a  digital  controller  can  not  be  used  because  of 
practical  limitations.  Therefore,  the  reduced  order 
feedback  controller  is  used  in  this  study.  It  is  assumed 
that  the  order  of  the  deterministic  ARMA  model  is  less 
than  or  equal  to  the  order  of  the  generalized  displace¬ 
ment  vector. 

M  <  m  , 

The  reduced  order  ARMA  model  (  2M<NM  )  faces  a 
potential  problem  since  unmodeled  high  frequency 
components  can  cause  a  frequency  aliasing  in  the 
parameter  estimation  and  control  law  design  pro¬ 
cedures.  Using  an  anti-aliasing  filter  is  one  possible 
approach  for  reducing  the  sensitivity  of  the  estimator. 

Since  the  contribution;  of  lightly  damped  modes 
dominates  the  dynamic  response  of  an  aeroservoelsuc 
system,  and  since  our  interest  is  focused  on  these 
lightly  damped  modes,  relatively  small  M  and  m  are 
selected  for  the  reduced  order  ARMA  model  and  the 
order  of  the  generalized  displacement  vector,  respec¬ 
tively.  In  this  study,  most  of  the  numerical  examples 
use  m  =  3  and  NG  =  8.  Therefore,  NM  becomes  30. 
However,  the  order  of  the  ARMA  model  2M  is  equal  to 
6.  Hence,  24  frequency  components  are  neglected  in 
the  adaptive  optimal  control  law  design  procedure. 

In  order  to  study  potential  problems  due  to 
unmodelled  dynamics  a  large  order  aeroservoelastic 
system  is  also  used.  A  tenth  order  sine  Butterworth 
low-pass  filter50  is  used  as  the  anti-aliasing  filter  in  the 
present  study.  The  number  of  structural  natural  fre¬ 
quencies  which  are  smaller  than  the  cutoff  frequency  of 
the  anti-aliasing  filter  is  assigned  to  the  value  of  M. 
For  this  example  with  the  unknown  random  external 
loads,  8,  16,  and  3  are  assigned  to  values  of  m,  NG,  and 
M,  respectively.  Thus,  eight  modes  are  included  in  the 
time  histories  computed  with  three  modes  in  the  previ¬ 
ous  case. 

Learning  Period.  A  short  learning  period, 
t  =  8 xM  xTc ,  is  used  for  obtaining  the  initial  estimates 
of  the  ARMA  coefficients  which  are  needed  for  initial¬ 
izing  the  Riccati  matrix  [PJ*  in  Eq.(34).  During  this 
period,  the  control  surface  is  excited  by  small- 
amplitude  white  noise.  At  the  end  of  the  learning 
period,  the  initial  condition  of  the  the  Riccati  matrix, 
[PJo,  is  calculated  from  the  discrete-time  algebraic  Ric¬ 
cati  equation,  Eq.(36),  using  Potter’s  method.51 

Caution  should  be  exercised  when  the  computer 
simulation  is  started  under  the  unstable  flight  condition 
because  wing  response  in  this  case  can  become  exces¬ 
sively  large  during  the  learning  period. 

Forgetting  Factor.  A  fundamental  property  of 
the  adaptive  controller  for  an  active  flutter  suppression 
system  is  its  ability  to  track  variations  in  the  aeroelastic 
system.  To  follow  changes  in  the  aeroelastic  system  it 
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.  is  necessary  to  discard  old  time  historyof  the  aeroelas- 
tic  response  information:  in  the  parameter  estimator.  It 
should  be  emphasized  that  tracking  a  system  with 
rapidly  varying  parameters  is  not  possible.  However,  a 
slow  time  varying  system  can  be  tracked  reasonably 
well.  It  is  assumed  in  this  study  that  aeroelastic  system 
parameters  change  slowly. 

Old  time  history  data  of  the  aeroelastic  response 
is  discarded  exponentially  in  this  study.40  When  the 
forgetting  factor  X*,  given  in  Eq.(33),  is  equal  to  one, 
all  data  have  the  same  weighting.  This  forgetting  fac¬ 
tor  is  used  for  the  time  invariant  aeroelastic  system. 
When  the  forgetting  factor  is  less  than  one,  recent  data 
are  given  more  weight  than  old  data.  Thus,  a  forget¬ 
ting  factor  of  less  than  one  is  used  for  the  lime  varying 
aeroelastic  system. 

From  this  brief  discussion  it  is  evident  that  the 
forgetting  factor  is  linked  to  the  rate  of  change  of  the 
system  parameters.  In  practical  implemetation,  changes 
in  the  properties  of  an  aeroelastic  system  can  be  associ¬ 
ated  with  changes  in  air  temperature,  pressure,  and 
speed  of  flight.  The  proper  forgetting  factor  for 
different  rates  of  parameter  changes  can  be  determined 
accordingly. 

Computational  Delay.  The  on-line  computation 
of  the  aeroelastic  system  parameter  estimation  and  con¬ 
trol  law  design  will  produce  a  computational  delay 
between  measurement  and  control  command.  In  this 
study  a  computational  delay  is  not  taken  into  account. 
It  should  be  noted  that  in  order  to  have  feasible  pracu- 
cal  implementation,  the  computational  delay  should  be 
less  than  or  equal  to  the  sampling  time  Tt. 

4.  The  Wing  Model  Used  for  Computer  Simulations 

The  structural  modelling  in  this  study  can  handle 
practical  wings  with  composite  skins.36  However  for 
the  computations  carried  out  in  this  study,  an 
aluminum/plastic  foam  wing  model  was  selected, 
because  the  results  of  the  computed  flutter  boundaries 
could  be  compared  with  experimental  results.  This 
cantilevered  rectangular  wind  tunnel  model  wing  with 
6%  circular-arc  cross  sections  and  an  aspect-ratio  of  5.0 
was  actually  built  and  tested  in  1959.52  It  has  an  alumi¬ 
num  insert  covered  with  flexible  plastic  foam.  Experi¬ 
mental  flutter  boundaries  and  natural  frequencies  for 
this  wing  are  available. 

Natural  frequencies  and  flutter  boundaries  for  this 
wing  model  were  computed  by  Guruswamy  et  al.  using 
XTRAN3S53.  and  XTRAN3S-AMES54  codes.  The 
planform  geometry  of  this  wing  has  been  modified  in 
the  flutter  suppression  studies  by  Guruswamy  et  al.4  A 
trailing  edge  control  surface  as  shown  in  Figure  1  was 
introduced  in  Ref.  4.  Five  translational  and  five  rota¬ 
tional  springs  are  used  in  this  study  to  physically  con¬ 
nect  the  wing  and  control  surfaces.  These  springs  are 
equally  spaced  along  the  hinge  line  between  the  wing 


and  trailing-edge  control  surfaces  as  shown  in  Figure  1. 
The  material  properties  of  the  model  wing  and  the 
additional  spring  constants  defined  in  this  study  are 
given  in  Table  1. 

Three  wing  rectangular  sections  as  shown  in  Fig¬ 
ure  1  are  used  for  the  LISSA  and  PCKFM  computer 
representation.  The  full  potential  CFD  code  developed 
by  Shankar  et  al.13  with  88x12x18  grid  is  used  to  obtain 
the  transonic  correction  matrices.  Steady  state  tran¬ 
sonic  pressure  distribution  on  the  wing  surface  for 
different  Mach  numbers  are  shown  in  Figure  2. 

5.  Results 

5.1  Flutter  and  Divergence  Boundaries 

Before  comparing  the  flutter  and  divergence 
boundaries,  natural  frequncies  of  the  cantilevered  rec¬ 
tangular  wing  obtained  from  the  on-line  parameter  esti¬ 
mation  technique  are  compared  with  the  experimental 
and  other  computational  results.  In  this  example  struc¬ 
tural  transient  time  response  was  calculated  using  the 
iterative  time  marching  structure/aerodynamic  integra¬ 
tion  scheme  with  the  assumption  (R  )4  =  {0}  and 
(Cw]  =  [0].  Natural  frequencies  from  the  experiment52  , 
a  finite  element  analysis53  ,  the  LISSA  code  with 
eigenanalysis,  and  the  LISSA  code  with  the  on-line 
parameter  estimation  technique  are  presented  in  Table 
2.  The  results  in  Table  2  are  calculated  for 
T,  =  0.001 656sec;  24  AR  coefficients  were  used  together 
with  45  sampling  points.  The  on-line  estimauon  of 
structural  frequencies  and  the  assumed  structural  modal 
damping  factors  for  the  cantilevered  rectangular  wing 
without  aerodynamic  loads  is  presented  in  Table  3. 
The  results  in  Table  3  are  calculated  for  the  same 
values  of  T,,  AR  coefficients  and  sampling  points  as 
used  in  Table  2.  It  was  also  assumed  that  the  modal 
damping  ratio  was  £  =  .0200,  i'=l,2,  •  ■  -  .5.  It  is  evident 
from  Tables  2  and  3  that  the  on-line  estimation  can 
yield  system  frequencies  and  values  of  modal  damping 
correctly.  In  Table  2  the  difference  bctw’een  the 
present  analysis  and  the  finite  element  analysis  of  Ref. 
53  is  that  mass  and  stiffness  properties  of  the  plastic 
foam  were  not  included  in  Ref.  53  analysis.  In  our 
analysis,  the  stiffness  and  mass  properties  of  the  flexi¬ 
ble  plastic  foam  are  included  in  the  calculation.  It  is 
evident  from  the  results  shown  in  Table  2  that  the 
stiffness  and  mass  properties  of  the  flexible  plastic 
foam  need  to  be  included  in  order  to  obtain  the  better 
agreement  with  experiment. 

Flutter  Boundaries.  In  the  estimation  of  the 
open-loop  flutter  boundaries  the  MA  coefficients  in  the 
ARMA  model  are  neglected.  Aeroelastic  stability  of 
the  wing  is  determined  from  the  eigenvalues  of  the 
estimated  matrix  (APJ  in  Eq.(29).  It  should  be  noted 
that  the  ARMA  model  in  Eq.(28)  is  based  on  the 
discrete  lime  system.  Aeroelastic  modal  dampings  and 
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*  damped  frequencies,  o;  and  a>dj,  for  the  continuous 
time  system  can  be  obtained  by  using  the  following 
equations. 

o/s-jjr  log,  (c/+-d/) 

1  -i  di 

**  =  r t8n  t 

Here,  c-t  and  dj  are  the  real  and  imaginary  part  of  the 
eigenvalues  obtained  from  the  estimated  matrix  (Ap). 
Flutter  boundaries,  for  the  aeroelastic  system  are  deter¬ 
mined  from  the  aeroelastic  modal  dampings  q,.  When 
Oj  <  0  the  system  is  stable  and  when  os  =  0,  it  is  on  the 
flutter  boundary.  When  o;-  >  0  the  system  is  unstable. 

The  effectiveness  of  the  parameter  estimation 
technique  based  on  the  ARMA  model  to  identify  the 
flutter  boundary  is  demonstrated  in  Table  4.  The 
results  in  Table  4  are  calculated  for  the  same  values  of 
Tt  and  AR  coefficient  as  used  in  Tables  2  and  3; 
qD  =  9.439  kPa  and  10  iterations  are  used  for  the  itera¬ 
tive  time  marching  structure/aerodynamic  integration 
scheme.  Flutter  boundaries  obtained  from 
eigenanalysis  of  the  aeroelastic  system  equation  in  first 
order  state  variable  form,  V-g  method  in  frequency 
domain,  and  the  on-line  parameter  estimation  technique 
applied  to  the  aeroelastic  transient  time  responses  are 
compared  in  Table  4.  In  this  table  subsonic  aero¬ 
dynamics  are  used  to  check  the  accuracy  of  the  itera¬ 
tive  time  marching  structure/aerodynamic  integration 
scheme.  The  agreement  between  these  three  sets  of 
results  is  quite  good.  Small  discrepancies  between  the 
eigenanalysis  and  on-line  parameter  estimation 
observed  in  Table  4  are  mainly  due  to  the  linear 
assumptions  for  fr?(/)}.fr(0),  and  {r(r)J  vectors  in 
time  interval  kTa<t<(k+l)Ta. 

Transonic  flutter  boundaries  for  the  cantilevered 
wing,  obtained  from  the  on-line  parameter  estimation 
technique  developed  in  this  study,  are  compared  with 
experimental  and  other  computational  results  in  Figures 
3  and  4.  Steady  state  pressure  distributions  for  various 
Mach  number  (Figure  2)  show  that  transonic  shock 
effects  start  to  appear  near  Mm  =  0.84.  Thus,  the  tran¬ 
sonic  approximation,  shown  in  Figures  3  and  4,  is  in 
reasonable  agreement  with  experimental  results  up  to 
Af.  =  0.84.  Between  Mach  number  0.84  to  0.90,  the 
shock  is  well  developed  on  the  wing  surface.  Some 
discrepancies  in  this  Mach  number  range  in  Figure  3 
and  4  may  be  due  to  the  mixed  subsonic  and  super¬ 
sonic  sub-regions.  When  the  Mach  number  exceeds 
0.90,  the  shock  will  reach  the  wing  trailing  edge  and 
supersonic  sub-region  will  become  much  wider  than  the 
subsonic  sub-region,  as  indicated  in  Figured.  There¬ 
fore,  the  subsonic  assumption  used  for  generating  the 
unsteady  transonic  aerodynamics  is  violated  and  it  can 
not  be  used  for  Mach  numbers  greater  than  0.90. 

Flutter  boundaries  with  the  steady  transonic  aero¬ 
dynamic  load  are  also  presented  in  Figures  3  and  4  to 


check  the  quasisteady  assumption  for  the  transonic 
correction  matrix.  It  is  evident  from  Figures  3  and  4 
that  the  quasisteady  assumption  for  the  transonic 
correction  matrix  gives  good  agreement  with  the  exper¬ 
imental  and  CFD  results  in  the  low  transonic  Mach 
number  range. 

Divergence  Boundaries.  A  transonic  divergence 
boundary  for  the  wing  can  be  obtained^  using  the  tran¬ 
sonic  aerodynamic  influence  coefficient  matrix  given  in 
Eq.(6).  From  Eqs.(2),  (5),  and  (6),  the  following  equa¬ 
tion  can  be  derived  to  include  only  structural  stiffness 
and  aerodynamic  stiffness  of  the  wing. 

( fro2]- [4,(0]  )fri(0)  =  {0] 

Divergence  boundaries  are  a  function  of  dynamic  pres¬ 
sure  qD,  and  it  can  be  obtained  from  the  determinant  of 
the  matrix  ( [co2J  -  (4,(r)] )  using  a  root  finder. 

Divergence  boundaries  for  the  cantilevered  rec¬ 
tangular  wing  using  subsonic  or  transonic  aerodynamics 
are  shown  in  Figure  5.  Comparing  Figures  2  and  5  a 
transonic  bucket  (  evident  also  in  the  flutter  boundary  ) 
occurs  at  a  lower  Mach  number  for  the  divergence 
boundaries.  In  the  present  case  the  transonic  bucket  in 
the  flutter  boundaries  could  be  directly  associated  with 
the  steady  transonic  aerodynamic  loads. 

5.2  Adaptive  Subsonic  Flutter  Suppression 

The  effectiveness  and  versatility  of  the  adaptive 
optimal  controller  is  studied  by  applying  it  in  two  flight 
regimes:  subsonic  and  transonic.  The  subsonic  flutter 
suppression  studies  are  used  to  learn  about  input 
parameters  for  the  digital  adaptive  optimal  controller. 

Random  External  Loads.  Robustness  of  the 
digital  adaptive  controller  was  tested  using  random 
loads.  The  computer  simulation  was  carried  out  using 
subsonic  aerodynamic  loads  obtained  from  assuming 
that  [T,  (»)]  =  (Tj(<»)]  =  [I],  The  random  loads  are  gen¬ 
erated  by  introducing  appropriate  random  changes  in 
the  angle  of  attack.  Maximum  change  in  the  angle  of 
attack  is  determined  by  using  a  suitable  combination  of 
gust  and  free  stream  velocities.  A  gust  velocity  of 
39.37  inch/sec  at  Mm  =  0,714  and  q0  =  9.763  kPa  pro¬ 
duces  approximately  one  inch  deflection  at  a  sensor 
position  indicated  in  Figure  1. 

The  generalized  equations  of  motion  are  based  on 
8  structural  modes  which  are  obtained  from  the  LISSA 
code  with  24  degrees  of  freedom.  The  aerodynamic 
influence  coefficients  evaluated  by  the  PCKFM  are 
based  on  the  same  number  of  degrees  of  freedom  used 
in  the  LISSA  model.  Aerodynamic  influence 
coefficient  matrices  are  tabulated  for  21  frequencies, 
where  the  highest  frequency  is  95.49 Hz.  Tabulated 
generalized  aerodynamic  influence  coefficient  matrices 
are  fitted  in  Laplace  domain  using  Roger’s  approxima- 
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rion  with  16  aerodynamic  lag  terms  which  are  equally 
spaced  between  0  and:95.49  Hz. 

A  sine  Butterworth  lowpass  filter  of  order  10  is 
used  as  an  anti-aliasing  filter.  Note  that  the  high  fre¬ 
quency  aliasing  effects  can  not  be  perfectly  removed 
with  the  anti-aliasing  filter.  It  should  be  emphasized 
that  the  parameter  estimator  (BUD  algorithm  combined 
with  covariance  resetting  technique)  which  we  use  is 
very  powerful  so  that  the  high  frequency  components 
can  be  identified  whether  the  anti-aliasing  filter  is  used: 
or  not.  When  the  anti-aliasing  filter  is  used,  estimated 
damping  in  the  high  frequency  components  are  much 
larger  than  for  low  frequency  components.  In  this 
research  one  anti-aliasing  filter  is  used  since  it  can 
remove  most  high  frequency  components  in  the  filtered 
aeroelastic  transient  time  responses,  in  an  adequate 
manner. 

The  time  step  of  Ta  =  0.000414  sec  is  used  in  the 
iterative  time  marching  struclure/acrodynamic  integra¬ 
tion  procedure.  For  the  parameter  estimator  the  sam¬ 
pling  time  is  taken  as  T,  =  0.00414  sec.  The  sampling 
time  T,  should  satisfy  the"  Shannon's  sampling 
theorem.55  Therefore,  the  Nyquist  frequency,  /mu, 
becomes  120.8  Hz,  and  the  first  three  aeroelastic  modes 
can  be  identified  using  the  parameter  estimator.  The 
ARMA  model  consisting  of  6  AR  and  6  MA 
coefficients  is  used  to  obtain  the  acroscrvoclasiic  state 
space  description.  Covariance  matrix  resetting  interval 
for  the  parameter  estimator  is  equal  to  12  sampling 
steps  (=  0.05  see). 

Weighting  factors  (Q)  and  rw  in  the  performance 
index  2,  given  in  Eq.(24),  are  assumed  as  [I]  and  30, 
respectively.  In  this  study  maximum  amplitude  of  the 
allowable  control  surface  deflection  angle  is  constrained 
to  be  less  than  4  degrees. 

The  following  computer  simulation  lasting  10 
seconds  is  made  of  four  important  time  periods.  The 
first  24  sampling  steps  constitute  the  learning  period. 
During  this  period  the  control  surface  is  activated  ran¬ 
domly  to  get  the  initial  estimation  of  the  aeroservoelas- 
lic  system  matrices.  A  maximum  amplitude  of  random 
control  surface  deflection  angle  is  limited  to  0.8°,  At 
the  end  of  the  learning  period  the  controller  is  engaged 
and  the  initial  condition  for  the  Riccati  matrix  is 
obtained  using  the  Potter’s  method.  In  the  second 
period,  from  0.1  sec  to  2.0  sec,  the  active  flutter  suppres¬ 
sion  system  controls  the  acceleration  at  the  sensor  posi¬ 
tion,  caused  by  the  random  control  surface  deflection 
angle  in  the  learning  period.  During  the  third  period, 
from  2.0  sec  to  8.0  sec,  random  variation  in  the  wing 
angle  of  attack  is  introduced.  The  robustness  of  the 
active  flutter  suppression  system  with  respect  to  the 
unknown  external  disturbance  can  be  verified  by  exam- 
ing  the  time  histories  of  the  controller  gain  vector. 
This  unknown  external  disturbance  is  removed  during 
the  fourth  time  period,  from  8.0  sec  to  10.0  sec.  A  con¬ 
stant  forgetting  factor  of  0.9999  for  the  parameter  esti¬ 


mator  is  used  throughout  the  computer  simulation  pro¬ 
cedure. 

Time  histories  of  acceleration  at  the  sensor  posi¬ 
tion  and  the  first  element  of  the  controller  gain  vector 
are  shown  in  Figures  6  and  7.  In  Figure  6  accelera¬ 
tions  at  the  beginning  of  the  second  period  are  rela¬ 
tively  large.  It  should  be  mentioned  that  estimated 
aeroservoclastic  parameters  at  the  beginning  of  the 
second  period  may  not  be  completely  converged  after 
24  time  steps.  It  is  evident  in  Figure  7  that  the  con¬ 
troller  gain  vector  is  robust  with  respect  to  the  unk¬ 
nown  external  disturbance.  During  the  third  period  the 
controller  gain  does  not  change  significantly,  and  when 
the  unknown  external  disturbance  is  removed,  the  con¬ 
troller  gain  converges  again  to  the  constant  value. 
Small  variations  in  controller  gain  can  be  observed 
before  and  after  the  third  period,  these  small  variations 
can  be  eliminated  if  the  longer  learning  period  is  used. 

Computation  times  required  for  the  parameter 
estimator  and  the  on-line  control  law  design  procedure 
are  shown  in  Table  5.  The  computational  delay  of 
0.002436  sec,  for  the  ARMA  model  of  order  6  on 
SUN3/280  computer  with  a  floating  point  accelerator  is 
smaller  than  the  sampling  time  T,  =  0.00414  sec.  There¬ 
fore,  the  digital  adaptive  optimal  controller  in  this  study 
is  a  feasible  for  practical  active  flutter  suppression. 

Time  Varying  Aerodynamic  Loads.  Time  vary¬ 
ing  subsonic  aerodynamics  are  obtained  from  the  table 
of  subsonic  aerodynamic  influence  coefficient  matrices 
as  follows: 

1)  Subsonic  aerodynamic  influence  coefficient 
matrices  arc  computed  at  different  Mach  numbers 
and  frequencies. 

2)  At  each  Mach  number  tabulated  subsonic  aero¬ 
dynamic  influence  coefficient  matrices  are  fitted 
with  respect  to  frequency,  using  Roger’s  approxi- 
mau'on. 

3)  Continuously  fitted  aerodynamic  influence 
coefficient  matrices  obtained  in  step  2)  are  fitted 
with  respect  to  Mach  number.  Between  two  adja¬ 
cent  Mach  numbers,  linear  interpolation  is  used  for 
simplicity. 

It  is  evident  from  Figure  6  that  high  frequency 
components  do  not  affect  the  performance  of  the  con¬ 
troller.  Therefore,  the  first  three  structural  modes 
without  the  anti-aliasing  filter  will  be  used  in  the  com¬ 
puter  simulations  which  follow.  Aerodynamic  influence 
coefficient  matrices  are  tabulated  at  11  frequencies 
between  0  and  47.75  Hz.  Twelve  Mach  numbers 
between  A/„  =  0.714  and  Af_  -  0,920  are  selected  to  gen¬ 
erate  the  time  varying  subsonic  aerodynamic  loads. 
The  Roger’s  approximation  with  8  aerodynamic  lag 
terms,  which  are  equally  spaced  between  0  and 
47.75  Hz,  are  used  to  fit  the  aerodynamic  influence 
coefficient  matrices  at  each  Mach  number. 
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A  computer  simulation  lasting  10  seconds, ‘which 
is  similar  to  that  described  in  the  previous  section  is 
shown  in  Figure  8.  The  main  difference  between  the 
previous  and  the  current  simulation  is  in  the  third 
period,  extending  from  2.0  sec  to  8.0  sec.  During  this 
time  period  the  free  stream  Mach  number  is  increased 
gradually  from  Af„  =  0.714  to  Af„  =  0.870,  and=therefore 
this  period  is  denoted  as  the  varying  parameter  time 
period.  During  this  period,  the  forgetting  factor  of 
0.9999  is  changed  to  0.85  to  follow  the  acroservoclastic 
system  changes.  The  flight  path  at  P„  =  20.21  kPa 
together  with  the  flutter  and  divergence  boundaries  are 
presented  in  Figure  9.  All  the  other  input  data  are 
identical  to  the  previous  computer  simulation. 

Time  histories  of  acceleration  at  the  sensor  posi¬ 
tion  and  the  first  element  of  the  controller  gain  vector 
are  shown  in  Figures  10  and  11.  In  Figure  11  con¬ 
troller  gain  follows  the  aeroservoelastic  system  changes 
successfully.  It  should  be  noted  that  the  robustness  of 
the  controller  gain  in  lime  varying  interval  depends  on 
the  forgetting  factor,  [Q],  and  rw. 

Free  stream  pressure,  P„,  for  the  second  flight 
path  in  Figure  9  is  assumed  to  be  27.36  kPa.  In  this 
case  free  stream  Mach  number  is  increased  gradually 
from  =  0.714  to  A/„  =  0.916.  Time  histories  of 
acceleration  at  the  sensor  position  and  the  first  element 
of  the  controller  gain  vector  are  shown  in  Figures  12 
and  13.  In  Figure  12  the  adaptive  flutter  suppression 
system  fails  around  time  t  =  4.5  see.  The  corresponding 
Mach  number  can  be  obtained  from  Figure  8,  and  is 
approximately  equal  to  =  0.8.  In  Figure  9  the 
second  flight  path  and  the  divergence  boundary  inter¬ 
sect  around  =  0.8.  This  is  another  example  of  the 
failure  of  acceleration  based  control  in  stabilizing 
divergence  type  instabilities.  To  suppress  both  the 
flutter  and  divergence  instabilities,  displacement  sensing 
is  used  along  the  second  flight  path,  and  results  are 
presented  in  Figures  14  and  15.  It  is  evident  that  both 
instabilities  are  now  suppressed.  Similar  results  have 
been  reported  by  a  number  of  authors.56,57 

5.3  Adaptive  Transonic  Flutter  Suppression 

Time  varying  transonic  aerodynamics  are 
obtained  from  tables  of  steady  transonic  and  unsteady 
subsonic  aerodynamic  influence  coefficient  matrices.  In 
addition  to  the  three  steps  described  in  section  5.2.2, 
the  tabulated  steady  transonic  aerodynamic  influence 
coefficient  matrices  are  also  fitted  linearly  between  two 
adjacent  input  Mach  numbers.  Thus,  the  time  varying 
transonic  aerodynamic  loads  can  be  obtained  with  rela¬ 
tive  ease. 

Transonic  flutter  boundaries,  transonic  divergence 
boundaries,  and  two  flight  paths  are  shown  in  Figure 
16.  Transonic  aerodynamic  influence  coefficient 
matrices  are  computed  at  A/„  =  0.714,  0.850,  0.873,  0.895, 
and  0.916.  The  same  input  data  given  in  section  5.2.2 
are  used  for  all  other  variables. 


Aeroservoelastic  response  with  acceleration  sens¬ 
ing  is  shown  in  Figures  17  and  18.  In  Figure  18  the 
controller  gain  is  increased  from  /  =  2.0  *«.  to 
t  =  6.5  sec  However,  from  t  =  6.5  sec  to  i  =  8.0  tec  the 
controller  gain  decreases.  In  Figure  8,  time  /  =  65  tec 
corresponds  to  A/..  =  0.87,  and  this  Mach  number  is 
approximately  the  center  of  the  transonic  bucket  in  Fig¬ 
ure  16.  It  is  evident  from  Figure  18  that  the  controller 
gain  is  increased  when  the  aircraft  penetrates  further 
into  the  flutter  region.  On  the  other  hand,  if  the  aircraft 
leaves  the  flutter  region,  then  controller  gain  will 
decrease. 

Aeroservoelastic  response  with  the  acceleration 
sensing  along  the  second  flight  path  is  shown  in  Figures 
19  and  20.  In  Figure  20  the  parameter  estimator  can 
not  follow  the  acroelstic  system  change  around 
/  =  5.5  sec.  Here,  time  r=  5.5  sec  corresponds  to 
A/..  =  0.83  in  Figure  8.  In  Figure  16  the  second  flight 
path  and  transonic  divergence  boundaries  intersect  each 
other  around  A/_  =  0.81.  Thus,  again  it  is  evident  that 
divergence  type  instabilities  can  not  be  suppressed 
using  the  acceleration  sensing.  The  aeroservoelastic 
response  with  displacement  sensing  along  the  second 
flight  path  are  presented  in  Figures  21  and  22.  The 
change  of  controller  gain  in  Figure  22  is  similar  to  that 
in  Figure  18.  Thus,  one  may  conclude  that  the  digital 
adaptive  optimal  controller  used  in  this  study  does  fol¬ 
low  system  changes  correctly. 

6.  Concluding  Remarks 

A  new  adaptive  optimal  control  methodology  for 
active  flutter  suppression  is  studied  for  a  number  of 
aeroclastic  systems.  This  controller  is  robust  with 
respect  to  the  unknown  external  disturbances  and  can 
be  applied  to  time  varying  flight  conditions.  The  same 
adaptive  optimal  control  algorithm  is  applied  in  both 
subsonic  and  transonic  flight  conditions. 

It  is  also  shown  in  this  study  that  the  divergence 
type  instability  can  not  be  adaptively  controlled  with 
the  acceleration  sensing.  Since  the  controller  perfor¬ 
mance  with  the  acceleration  sensing  strongly  depends 
on  the  relative  positions  of  flutter  and  divergence  boun¬ 
daries,  it  is  important  to  design  the  controller  and  the 
structure  simultaneously  in  an  integrated  manner  so  as 
to  avoid  such  difficulty. 

The  ARMA  model  with  a  parameter  estimation 
technique  is  successfully  implemented  to  yield  tran¬ 
sonic  flutter  boundaries  of  the  wing  structure.  This 
methodology  yields  flutter  boundaries  very  efficiently  in 
time  domain,  and  thus  only  a  small  number  of  aeroelas- 
tic  transient  time  responses  are  needed  in  order  to 
obtain  reasonable  results. 

A  simple  methodology  for  obtaining  three- 
dimensional  unsteady  transonic  aerodynamics  in  the 
time  domain  is  presented.  The  transonic  bucket  is  suc¬ 
cessfully  reproduced  and  predicted  by  this  simple 
approximate  method.  Since  the  unsteady  computations 
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are  based  on  unsteady  subsonic  aerodynamics,  transonic 
acroclastic  responses  in  the  time  domain  can  be 
obtained  at  low  computational  cost.  The  most  time 
consuming  ingredient  in  this  approximate  method  is  the 
computation  of  transonic  aerodynamic  influence 
coefficient  matrices. 
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Table  1.  Material  Properties  of  the  Cantilevered  Rec¬ 
tangular  Wing  Model 


Material  Prooertiei  1 

Young’s  Modulus 

9847878.6  psi 

Poisson's  Ratio 

03528521 

Density 

0.1116599  lbmlin' 

Translational  Spring  Constant 

ljt9  Ibf  tin 

Rotational  Sorine  Constant 

l.e  1  Ibf  -intrad 

Table  2.  Natural  Frequencies.  (Hz)  of  the  Cantilevered 
Rectangular  Wing 


Mode 

Exceriment* 

FEM*» 

LISSA 

LISSA* 

LISSA*i 

1-Sl 

14.29 

13.21 

15.18 

15.18 

15.20 

2-nd 

80.41 

67.32 

78.04 

78.04 

75.94 

3-rd 

89.80 

76.97 

94.46 

94.46 

9658 

4-th 

— 

203.6 

250.5 

250.5 

225.8 

5-th 

— 

203.9 

266.0 

266.0 

274  7 

•  Ref.52;  ••  Ref.53;  t  On-line  parameter  estimation  technique  u 
applied  to  the  free  vibration  of  the  wing:  *  Cantilevered  wing 
with  control  surface. 


Table  3.  Natural  Frequencies  (Hz)  and  Modal  Damping 
Factors  of  the  Cantilevered  Rectangular  Wing 
with  Control  Surface 


Mode 

Eigenf 

Estimated  CO.  * 

Estimated  C  • 

1-Sl 

15.20 

15.20 

.0200 

2-nd 

75.94 

75.94 

.0200 

3-rd 

96.58 

96.58 

.0200 

4-th 

225.8 

225.8 

.0200 

5-th 

274.7 

274.7 

.0200 

■f  Obtained  from  the  LISSA  code  using  eigenanalysu',  •  Calcu- 
CO/;  _  Oi 

ljicd  using  CO.  =  —  ~  *nd  C;  =  — .  Hoc,  CD*,  ukJ  O, 

*  a, 

are  estimated  from  the  on-line  parameter  estimation  with  the  AR 
model. 


Fig.l  Cantilevered  Rectangular  Wing  with 
Trailing-Edge  Control  Surface 
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Tabic  4.  Flutter  Boundaries  of  the  Cantilevered  Rec¬ 
tangular  Wing  at  Mm  =  0.714  from  the 
Eigenanalysis,  V-g  Method,  and  On-Line 
Parameter  Estimation  Technique 


Table  5.  Execution  Time  (seconds)  for  the  On-Line 
Computations  on  SUN'3/280  Computer  with 
Roaring  Point  Accelerator 


Order  of 
the  ARM  A 
Model 


lie?  nive 

Rieciti 

Solver* 


Compuu- 
tional 
Deli 


Distribution  on  the  Wing  Surface 


Flutter  Frequency  (Hz)  Dynamjc  Pressure  (Pa) 


:  Approximation  *  ■  :  Experiment 


Acctlrrarinn 


Arctlrnnon  ( inrhiifc * 


-  : Flutter  Boundaries  -  :  Path  2 


Mach  Number 

Fig.9  Computer  Simulation  Paths  for  Subsonic  Acroscrvoclasticity 


rs  ix  ijb  i«  **  Me  in  lie  •*  *  •*  in  t  is  ix  c  t  •«  « x  c  k  ix  ic  i«  ik 

Time  (  *10  *ec  )  Time  (  *10  yt  1 


Mg-10  Subsonic  Acroservociastic  Response  along 
Path  1 


Fig.1 1  Controller  Gain  for  Subsonic  Atroscrvo:Ia>u 
city  along  Path  1 


Acceleration 


Acceleration*  inch  arc* 


Acceleration 


